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Chapter 1 


Introduction 


Arguably every thesis concerning superconductivity starts with a few words about the great 
discovery by Heike Kamerlingh-Onnes, and this one shall not be different. The measurements 
of 8th April 1911 left an imprint in the world of physics, remarkable not only for the huge 
impact that the zero-resistance flow of charge has on the modern world, but also, or maybe 
even more for the gargantuan struggle of theoretical physicists aiming to erect a microscopic 
theory describing this baffling phenomenon. The epochal Kwik nagenoeg nul [1] reported by 
Kamerlingh-Onnes raised a puzzle with which giants like Bloch, Bohr, Heisenberg or Born 
wrestled and failed. Even Feynman, albeit concluding that the correct theoretical description 
of superconductivity is unachievable using perturbation theory, ascribed his low amount of 
published work in the 1950’s to his endeavor in finding the correct approach |2]. It took 46 
years and efforts of John Bardeen, Leon Cooper and Robert Schrieffer to finally arrive at 
the correct microscopic description of superconducting state, termed since as BCS theory [3]. 
Another big leap in understanding exotic phases of matter, which happens to be connected 
to the present thesis, was ignited by the discovery of the quantum Hall effect by Klaus von 
Klitzing [4]. As a result, the celebrated Landau paradigm of classifying phases of matter by 
symmetry breaking and order parameters became inadequate. Across the superconducting 
transition, it is the gauge symmetry that is said to be broken (although it is somewhat of a 
loose statement [5]), and an order parameter acquires a non-zero value. There is no local order 
parameter that can describe the quantum Hall state, nor is there a spontaneously broken 


symmetry when transitioning to such a state. We can rather say that this state is topologically 


distinct from other phases of matter, in the sense that some of its properties (e.g. number of 
gapless states dispersing on the edge of the sample) are invulnerable to smooth perturbations 
of the system like irregularities or disorder, as long as the single particle gap is not closed. 
It is the twisting and knotting of electron wavefunctions in momentum space |6] that is 
responsible for the existence of celebrated gapless edge states, when topologically disparate 
materials are connected — a manifestation of bulk-boundary corespondence. Whether they 
are end states in the Su-Schrieffer-Heeger model [7], helical edge states of the quantum 
spin Hall effect [8,9], or chiral Majorana modes on the edge of a semiconducting film in a 
two dimensional heterostructure [10], they have received an enormous amount of attention 
recently, especially since 2016, when John Michael Kosterlitz, David Thouless and Duncan 
Haldane were awarded the Nobel Prize in Physics for theoretical discoveries of topological 
phase transitions and topological phases of matter. Both superconductivity and topology 
are important in the present thesis, although it is the presence of impurities and bound 
states associated with them, that permiate throughout this work. On the one hand, there 
are scalar impurities, i.e. those which do not break time reversal symmetry and hence do 
not lift the Kramers degeneracy, which were shown by Anderson [11] to not alter the critical 
temperature of a superconductor. However that may be changed, when a disordered system is 
strongly correlated [12]. Additionally, non-magnetic impurities will not host any bound states, 
although this statement is not true for superconductors with pairing symmetry different 
than s-wave [13] and in the case of quantum dots. On the other hand, magnetic adsorbates 
locally break Cooper pairs and induce a pair of bound states in the superconducting energy 
gap. The focus on the smallest possible additives and their local influence on the host would 
not make sense without an appropriate experimental technique. The great invention by 
Gerd Binnig and Heinrich Rorher |14| (who shared the 1986 Nobel Prize in Physics with 
Ernst Ruska) revolutionized the domain of surface imaging, and provided an exceptional 
device — Scanning Tunneling Microscope (STM), which to this day provides astonishing 
images of surfaces and enables the study of various phenomena manifesting in the local 
density of states. The preceding phenomena appear throughout this work, and are tied 
together by one more - the spin-orbit interaction. The first contact with it is usually the 


introductory course on quantum mechanics, where one learns that it is a relativistic effect, 


which results in splitting of the atomic energy levels. In condensed matter physics, we are 
more used to the idea of electronic bands. Nevertheless it is convenient to consider an electron 
moving in the electric field stemming from the crystalline potential. In its rest frame, the 
electron feels an effective magnetic field, proportional to its momentum, and it is this effect 
that leads to energy splitting of the spin sub-bands. The spin-orbit interaction is odd in 
momentum, and the first to show that it arises in systems without an inversion symmetry was 
Dresselhaus [15] and Bychkov and Rashba [16]. To visualize the physical implications of the 
e.g. Rashba effect, one can formulate an analogy to the Magnus effect. A spinning object's 
trajectory will be deflected due to a force coming from a pressure difference. Analogously, the 
electrons moving through a solid may also be deflected, depending on their intrinsic angular 
momentum, resulting in spin Hall effect. Recently, spin-orbit coupling in solids receives much 
attention, not only because it is responsible for a surge in development of a new field dubbed 
spin-orbitronics, but also due to the emergence of many novel phases of matter which require 
the presence of this interaction [17]. This thesis focuses on bound states induced by the 
presence of magnetic impurities in superconductors and with the influence that various types 
of spin-orbit interactions might have on them. In the main part the theoretical background 
concerning most important phenomena will be laid out. Next we present a discussion of 
methods used to examine the selected physical systems. A summary of the findings from 
the papers co-authored by the author and closely related to the present thesis will serve as a 


conclusion. 
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Chapter 2 


Theoretical background 


2.1 Fundamentals of superconductivity 


As mentioned above and realized by many physicists in the 20th century, the correct description 
of the superconducting state cannot be approached by perturbative methods. It was Leon 
Cooper, who explicitly showed that it is the case |18|. Assuming that two electrons forming 
a singlet state interact above a filled Fermi sphere, Cooper showed that the energy (relative 


to the sum of energies at the Fermi sphere) of such a pair will become 
ui 
E€ = —2wpe V9»). (2.1) 


We note that this energy is negative, hence this state is a bound state. The Debye frequency 
wp refers to the idea of Fröhlich, that the attractive interaction between electrons might 
originate from electron-phonon interactions |19]. Building on these facts BCS proposed a 


variational wavefunction, which is a coherent superposition of Cooper pairs: 


IV scs) = | | (ua oe )10), (2.2) 
k 


with zero being the vacuum and uz, Uk some complex coefficients, on which we will elaborate 
below. One can see that there are many peculiar properties of this wavefunction, like the fact 


that it is a superposition of states with different total number of electrons. Knowing that 
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only electrons with energies below hwp are important in the problem, and that the Cooper 


instability arises from the scattering of electrons with opposite spin, we can write 


H = ` Erc! Cho + Dp Vie CIC jay Ce, Chet (2.3) 

ko k,k' 
and treat Vi, as an effective, attractive interaction (Vi; = V < 0) between electrons of 
interest (with |£&| < hwp). We will also focus on the isotropic (or s-wave) gap and write 
A, = A. Introducing the mean-field approximation, justified by the coherence length € in 
BCS superconductors being of the order of 1000A, we argue that the difference between an 


operator and its expectation value is small 
(A — (A))(B — (B)) ^ 0, (2.4) 


so we can say that 


AB c (AB + A(B) — (A)(B). (2.5) 


We then identify A — CIN UP and B = c jcy4, to finally express the reduced BCS 
Hamiltonian (with the constant term (A) (B) omitted): 


H zs ` Erc! Cho + A ` CIE e + A* ` C-k} Ck; (2.6) 
k 


ko k 


where 
A = V (ciet) 
(2.7) 
A= AC a) 
The non-zero value of A informs us about the superconducting transition, and as we will 
shortly see, its magnitude reflects the energy gap in the electronic spectrum. To analyze 


the Hamiltonian (2.6) we will employ the Bogoliubov-Valatin transformation, which can be 


summarized in the following identity: 
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in which we introduce new fermionic operators 2 which annihilate (create) Bogoliubov 


quasiparticles — Bogoliubons. To state the matter in a more picturesque fashion — an electron 
is represented as a superposition of a quasiparticle and a quasihole. We demand that the 
BCS Hamiltonian Eq. (2.6) is diagonal in the space of those operators. The explicit form 
of the operator space rotation matrix (2.8) can vary in literature. Common is the fact that 


when plugged into the BCS Hamiltonian it will yield a familiar formula: 


A 


UkUk = TS >> 
"EX 


which when combined with the result of enforcing fermionic anitcommutation relations of 


(2.9) 


Operators Yk 


{ot vet = 1, (2.10) 


gives the following relations: 


(2.11) 


Having those identities, we can insert them into the Hamiltonian expressed with the Bogoliubov 


quasiparticle operators, to see that it acquires a simple form: 


H = M Ey yee: (2.12) 

ko 
where the energy in the diagonal basis Ej = Je? + A2. We can see from (2.11) that in the 
normal state (A — 0), the Bogoliubons represent holes for energies less than the Fermi level, 
and particles for £& > Ep. The onset of superconductivity, embodied by the non-zero value 
of A, leads to a gap in the spectrum — the energy of an excited quasiparticle has to be at 


least A. We will now examine the gap equation 


A Ex 
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which can be obtained via insertion of the Bogoliubov-Valatin tranformation into Eq. (2.7), 


and study its limiting behavior. First, if we let T — 0 we are left with the following equation 


wp 
de 
1=Vp(0 —————. 2.14 
at T VFA pg 
0 
After careful evaluation, the above leads to a familiar result 
A = 2wp eV» (2.15) 


in agreement with the Schródinger equation approach used by Cooper (cf. Eq (2.1)). On the 
opposite side, when we send the gap to zero A — 0, we can find the critical temperature Tc. 


Again carefully evaluating the gap equation we arrive at 


To = — — ev» 2.16 
C x^ kp € ’ ( ) 
l (T = 0) 
where y is the Euler constant. Evaluation of LT leads to the most celebrated result 
BİC 
of the BCS theory, stating that 

A(T — 0) 
x~ 1.76 2.17 
ETE (2.17) 


is universal for conventional superconductors. 


Yu-Shiba-Rusinov states 


As every material is fundamentally impure, the study of impurities in both metals and 
superconductors is natural. When considering magnetic impurities in metals, one can dwell 
on whether the Kondo effect [20], or the Anderson impurity model [21] is the apex of 
theoretical efforts. The first studies examining a spherical impurity in superconducting hosts 
were reported by Fetter [22]. It was quickly realized that a non-magnetic impurity does 
not give rise to bound states with subgap energy, or rather, they do exist only in a window 
of about 1073A from the gap edge, thus they are irrelevant [23]. The above statement is 
modified in the context of 'artificial atoms’ - quantum dots (QD), which can support in-gap 


bound states when proximitized to a superconducting lead. When the charging energy of 
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the dot is smaller than A, there is a bound state composed out of a superposition of empty 
and doubly occupied state [24]. In the opposite regime (more common), i.e. U >> |A|, if 
the dot hosts an odd number of electrons, a single quasiparticle bound to the spin on the 
QD emerges. This state is usually called Yu-Shiba-Rusinov state, in analogy to the case of 
a classical magnetic impurity in a superconductor |25]. One can see that in the context of 
QDs, there is not much difference between the two regimes in terms of spectroscopy. This 
difference is however rather large when considering impurities, as we will now find out. The 
late 1960s were the time when the topic of bound states in s-wave superconductors induced 
by magnetic impurities was first considered by Luh Yu [26], Hiroyuki Shiba [27] and Anatol 
Rusinov [28]. The main caveat is that the impurities were assumed to be classical spins. It 
is probably put in the most elegant way in Shiba's paper, where he claims that having a 


Hamiltonian for the interaction of conduction electrons with a localized spin: 
H = — ` d ocw: S, (2.18) 


one can assume that the strength of the interaction J — 0 and the impurity spin S > oo, 
keeping the product JS finite. In this way we omit any quantum mechanical dynamics of spin. 
Despite the three aforementioned authors employing different approaches, the conclusion is the 
same: if there is a time reversal breaking, localized potential in contact with a superconductor, 
there will be a pair of single particle bound states inside the superconducting gap. First 
experimental results confirming this statement were reported in 1997 when the IBM Almaden 
group examined various magnetic and non-magnetic adsorbates on the surface of Niobium [29]. 
The STM spectra collected around the magnetic impurities revealed two bound states at 
energies below the BCS gap, which were asymmetric in spectral weight — a consequence 
of particle-hole symmetry breaking Coulomb potential and imbalances in the normal state 
conductance [30]. We will now embark on finding the energy of such bound states, starting 


with the familiar BCS Hamiltonian, or rather, the Hamiltonian density 


H = £y74,09 + A409 — JSÓ(T)T90;, (2.19) 
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in the Nambu basis V — (Vis, s Vrs =el u)T. In (2.19) the tensor product between 
matrices 7 and c is implied. These are the Pauli matrices which operate in particle-hole and 
spin space respectively. To make things simpler we can now split the Hamiltonian into two 


separate ones for spin-up and spin-down (while dropping the identity matrix To) 


H4 = &yT, + Ar, F JSÓ(r). (2.20) 


We now write the Schródinger equation with the impurity term isolated on the right hand 
side: 


[E — &r. — Ar]i(r) = JS9(r)u(0), (2.21) 


and proceed by employing a nifty trick. We first make use of the Fourier transform Yk = 


fA Ons e‘kru(r) to get rid of the Dirac delta, which gives: 
[E — €x, — Ar] = FISY(0). (2.22) 


Now we multiply from the left by [E — £y; — At,]~' and use the inverse Fourier transform 


to obtain an equation for the spinor at the position of the impurity, which is the origin: 


vOs f am s lE + bers + Are], (2.23) 


— oo 


where we have explicitly written the inverse matrix, an easy task considering we are working 


with 2x2 matrices. Looking at (2.23) we observe that we have two integrals to solve: 


T AÈ 1 
T dk A 
T = i (228 E -g EA. (2.25) 


We first make use of a usual substitution fA ons = = vo f d£y, where vo is the density of states 


at the Fermi energy. Now the integrand in Eq. (2.25) is odd in €% so the integral vanishes. 
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The integral (2.24) can be solved using contour integration, seeing that we are interested 
TVO 


V 


in subgap energies E < A, and we obtain 7; = — . Using this result we can now 


express equation (2.23) as: 


T E v(0) — 0. (2.26) 
Solving for energy we obtain the formula for the energy of Yu-Shiba-Rusinov states 


1-a? 
Eysn = D EEE (2.27) 


where in both of the above equations we have introduced the coupling constant as a = JS. 
By examining Eq. (2.27) we can see that the energy of a Shiba state will always be smaller 
than the the superconducting gap. Additionally, a special value of the strength of impurity- 
substrate coupling can lead to the energy of YSR states being zero. This is the critical 
coupling, at which the ground state of the system changes — a manifestation of the quantum 


phase transition, which we will now study in more detail. 


Quantum phase transition 


Since in general we focus on classical impurities with S >> 1, complete screening of this 
moment is not possible, hence the implications of the Kondo model, albeit admittedly 
fascinating [31,32], can be discarded. We will therefore follow the reasoning of Balatsky et 
al. [13] and references therein. In the weak coupling regime, i.e. when the coupling constant 
J is smaller than the critical value Jc = 1/795, the ground state is the true BCS state 
in the presence of the impurity potential (cf. Eq. (2.2)). If we examine the net spin of 
conduction electrons in the ground state we will naturally see that it is null (Wo|$_7|Wo)=0. 
The first excited state in this regime is the quasiparticle with energy Eysg and we label it 
as |V.,) = 4! |V9), where ^ is the quasiparticle creation operator that we have encountered 
in section 2.1. It is a fermionic quasiparticle, thus it is now not surprising that the net 
spin is (V,,[SZ|V.,) = —1/2. When the impurity coupling reaches the critical value, Ey sg 
becomes a zero energy state and the ground state becomes unstable. Above Jc it is more 
energetically favorable to induce the unpaired spin rather than to form a pair. Thence in the 


strong coupling situation J > Jc, |Wo) and |V.,) exchange roles, and now the ground state of 
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the system (with net spin —1/2) is a non-BCS state, because of the presence of an unpaired 
single-particle state. The additional spin-down quasiparticle (assuming antiferromagnetic 
alignment and spin-up impurity) locally supresses the pair potential A, even resulting in its 
sign change in the strong coupling regime. This is reminiscent of the 0 — 7 shift in Josephson 
junctions, manifested by the reversal of the supercurrent [33], and indeed the sign reversal 
of the order parameter in the vicinity of the impurity is sometimes dubbed ‘r phase shift’. 
While the case of Josephson current reversal is well established and understood, origins of 
the sign change of A(r — rg) remain mysterious, although numerical simulations suggest it is 
rather general. An analytic estimate was put forward in a recent work [34], which predicts 


the change in the gap function at the impurity site for J = Jc to be 


A 


ONE = Timp) ~ ON 


(2.28) 
It is also worth noting that when treated self-consistently, the supression of A around the 
impurity contributes to the free energy, and the critical value of the impurity coupling is 
shifted downwards. Naturally a question arises whether the presented reasoning reflects itself 
in any form in empirical results. Looking at the expression for Eysg Eq. (2.27), only the 
coupling J is a parameter that can be experimentally manipulated. The most common way 
is to depend on the discrete changes in J determined by the adsorption sites [31,35]. A new 
approach is to use the forces between the STM tip and the impurity (approximated by the 
Lennard-Jones potential) and continuously control the coupling, as presented for the first 


time in reference [34]. 


2.2 ‘Topological matter 


Topology as a branch of mathematics, which studies how shapes can be smoothly deformed into 
each other, was adapted in condensed matter to describe topologically equivalent Hamiltonians. 
It turns out that just as geometric objects, which are topologically equivalent have the same 
topological invariant, translationally invariant Hamiltonians, connected by adiabatic changes 


preserving the energy gap, can be called topologically equivalent. This naturally leads to 
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the fact that if the excitation gap closes upon such deformation, the topologically distinct 
systems are separated by a phase transition. If we are concerned with Hamiltonians labeled 
by the crystal momentum, the information about invariants must be contained in their 
wavefunctions. It was sir Michael Berry, who found that aside from the familiar dynamical 
phase, eigenstates of a Hamiltonian depending on slowly varying parameters pick up an 
additional, geometrical phase, now widely recognized as the Berry phase [36]. Exploiting the 
ambiguity in the solutions of Hy|u,(k)) = E,(k)|u,(k)), we can write 


[us (K)) — ei? 9 u, (K)). (2.29) 


The rate of change of the wave function in k space is defined as 


An(k) = (Qus (K)| V k]us(K)), (2.30) 


and bears the name Berry connection. It is an analogue of the electromagnetic vector potential 


and by virtue of the transformation Eq. (2.29), it also changes in a familiar way 


Anlk) —. An(k) — Vkos (Kk). (2.31) 


A gauge invariant quantity that we can construct is the Berry curvature, again in analogy to 
electromagnetism JF, (k) = Vp x A,(k). It turns out that the integral of the Berry curvature 


over the two dimensional Brillouin zone is the Chern number 


1 
C, = — f dk, dky Fnk), (2.32) 
BZ 


2m 
which is one the most fundamental integer invariants. In the above equation the superscript 1 
reflects the fact that it is the first Chern number, closely related to the quantum Hall effect, 
2 
e 
where the conductance is given by oz, = n—, with integer n being the total Chern number 


h 
of the occupied bands 


n=C= y C}. (2.33) 
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It was shown by Simon [37], that the language of Berry phase and the Chern number is 
equivalent to the famous TKNN invariant, formulated by Thouless, Kohmoto, Nightingale 
and den Nijs in the seminal 1982 paper [38]. There are different invariants for different classes 
of symmetries. Depending on the presence or absence of time reversal, particle-hole (charge 
conjugation) and chiral symmetries, topological states fall into one of the 10 symmetry classes, 
summarized in a "periodic table’ [39]. One additional remark is due as we segue to the 
description of the quantum spin Hall effect in the next paragraph. The total Berry curvature 


is odd under time reversal 


>, Falk) + — », Falk), (2.34) 

En<EF En<EF 
therefore the total Chern number of a system, which preserves time reversal symmetry is 
C = 0. It is no surprise, as the quantum Hall effect requires an externally applied magnetic 
field. However, the quantum spin Hall phase is an example of a topologically non trivial 
phase utilizing spin-orbit interactions, which do not break the time reversal symmetry. We 


will now shortly describe it. 


Quantum spin Hall effect 


Before introducing the ideas from the seminal paper by Kane and Mele [9], it is useful to 
mention the general approach taken by Haldane to describe a two dimensional topological 
insulator in a graphene lattice [40]. The continuum description of graphene, with the 


Hamiltonian linearized around the K or K’ point is 
H = —Hhwp(e37,0; + 0,0,), (2.35) 


with up = m where a is the lattice constant and t nearest-neighbor hopping. Pauli matrices 
c and T operate in the sublattice and valley (K and K’ point) spaces respectively. We have 
also already expressed the momentum measured relative to the Dirac points as q > —iA V. 
The celebrated linear spectrum follows from Eq (2.35) E(q) = hvrp|q|. At the Dirac points 


(q — 0) the inversion and time reversal 7 symmetries protect the degeneracy. The gap 
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can be opened by breaking those symmetries, e.g. by introducing a mass term proportional 
to oz, which would describe a situation with two different atoms in the unit cell, but this 
would give us a trivial insulating state. Haldane imagined introducing a second-neighbor 
hopping term, with different signs depending on the direction of hopping from the first to 
second neighbor. This would stem from some staggered magnetic flux, which breaks time 
reversal symmetry, but is equal to zero through the unit cell. The continuum version of the 
Haldane term acquires the form 


Hey = MAOzTz, (2.36) 


with opposite sign of the mass term at K and K’. This gapped system is a topological 
insulator, with the Hall conductance te?/h, depending on the sign of the Chern number, 
which can be +1. Solved in a semi infinite ribbon geometry, Haldane model will show a 
gapless state dispersing through the gap. This is the chiral edge state, which may cross 
the Fermi level multiple times with different group velocities. What remains robust is the 
difference between the right and left moving modes, which will be always equal to C = +1. 
The idea of Kane and Mele was to construct two copies of the Haldane model — one for each 
spin direction s, = +1. With addition of spin, a time reversal invariant term emerges, which 


couples spin and orbital degrees of freedom 


HKM Tz ASOOzTzS;. (2.37) 


In the lattice version of the Kane-Mele model there appears again a second neighbor hopping 
term, but it additionaly distinguishes between two spin species, for which this term has 
opposite signs. Because it is a doubled Haldane model, there are also two modes dispersing 
through the gap, however this time these are counterpropagating helical states with the spin 
locked to the momentum direction. Additionally, those states are immune to localization 
and will travel through a disordered medium without backscattering. One can show that 
the off-diagonal elements of the scattering matrix will vanish, as long as the time reversal 
symmetry is preserved. The problem remains, that the Chern number will be equal to zero 
in systems with time reversal invariance. If however the spin projection is a good quantum 


number, the difference between spin-up and spin-down Chern numbers will be C4 — C, = £2, 
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and can serve as an invariant. In reality the spins will be mixed and this recipe is no longer 
valid. Kane and Mele introduced a different topological index - the Z» invariant — in another 
seminal paper from 2005 [41]. The introduction of this Z number is deeply involved, and 
we will retain to a comment about the invariant in centrosymmetric systems. If there is an 
inversion symmetry, a Kramers pair of eigenstates at time reversal invariant momenta T; 


share an eigenvalue £; of the parity operator P 


Plu; (Ti) = élu (T), 
Pug (F2) = &lu; (P2) 


(2.38) 


with £; = +1 and the roman numerals label the Kramers partners. We can then evaluate the 


Zə index by computing 


ee IE (2.39) 


which probes if the edge states dispersing through the gap connect pairwise at the time 
reversal invariant momenta. If not, the invariant is (—1)" = —1 and the system is in the 


quantum spin Hall phase. 


Topological superconductivity 


Topological arguments are naturally extended from insulators to superconductors. After all, 
there is a gap in the spectrum of a superconductor, and states below the Fermi level are 
occupied. One can therefore define topological numbers for the occupied bands. Topological 
superconductivity is defined as a state with a fully opened gap, a non-zero topological invariant 
and the absence of gapless bulk excitations. Just as in the case of topological insulators, 
there will be gapless states only on the boundary between such a system, and a topologically 
trivial one. Depending on the bulk topology, the dispersion of those gapless states may be 
different; in the following paragraphs we wil focus on a nodal state with flat dispersion, and a 
chiral one, with a single state dispersing through the gap, similar to the quantum Hall state. 
The quasiparticle excitations in superconductors are a mixture or electrons and holes, so an 


eigenstate of the Bogoliubov de-Gennes Hamiltonian with energy E, has a partner with energy 
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—E. We see that a quasiparticle is essentialy the same as the antiparticle. A zero-energy 
fermionic excitation which comes to mind in this context is the Majorana fermion. Predicted 
in 1937 by Ettore Majorana, as a real solution of the Dirac equation, with a particle being 
its own antiparticle — a self-conjugate Dirac fermion [42]. In condensed matter, there are no 
real Majorana fermions known from high energy physics, rather there are emergent collective 
excitations of electrons [43]. In topological superconductors, the edge excitations are gapless, 
and their low energy description satisfies the massless Dirac equation, therefore the zero 
energy edge states are Majorana quasiparticles. The tremendous activity in researching the 
possible routes of obtaining a solid level of control of Majorana zero modes, comes from the 
prediction that they would serve as nonlocal qubits, because of their non-abelian statistics. 
The number of review papers elaborating on various means of inducing those exotic states in 
condensed matter systems is a testament to their relevance |44-48]. In the following we will 


shortly describe two types of systems, relevant to the contributions presented in Chapter 4. 


Chiral superconductivity 


The first demonstration of a chiral p-wave superconductor was put forward by Read and 
Green, with the gap structure explicitly breaking time reversal symmetry Ao(k, + iky) [49]. 
This state supports Majorana zero modes on the edges and in the vortex cores. A simpler 
realization was desired, and it turned out that Dirac fermions with s-wave pairing will also 
produce the desired outcome [50]. It was challenging at the time to find a condensed matter 
system realizing this proposal, but when it turned out that a surface state of a topological 
insulator is a Dirac fermion, Fu and Kane proposed, in the groundbreaking paper, to put 
such a surface state in proximity to s-wave superconductor [51]. The Hamiltonian of this 


system becomes 
H= (2.40) 
and hosts Majorana zero modes in vortices and on the edges. Another breakthrough 


demonstrates the important role of the Rashba spin-orbit interaction [10,52]. In the basis 


of We = (Cet Ck| c kt "m pe the Hamiltonian of a system with s-wave pairing, Zeeman field 
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and Rashba spin-orbit interaction is the following 


Ek - gk: € — ip Hio; iAc 
Hp = £ (2.41) 
—iAoc, —Ek + gk: C" + ug Hio, 
with gk = 2A(k,, —k,), present when the z + —z symmetry is broken. How this state is 
similar to a chiral p-wave superonductor can be seen by performing a unitary transformation 


HP = DH, DI, with 


D = — A (2.42) 
The dual Hamiltonian is now 


A-ppHzo,  —iEkOy — igk ` oo 
HP = HBHZ kOy — (Gk y (2.43) 
lEkOy + 1gk ` 90, —A + ppHzo, 
and we see that the spin-orbit coupling now serves as a triplet pairing component, similar to 


the Read & Green’s proposal [49]. After a critical value of magnetic field is applied 


upH, > V/e(0) + A3, (2.44) 


a topological gap opens and an edge state with linear dispersion E ~ ck traverses the gap. 
Applying a strong field can of course destroy s-wave superconductivity, but the difficulty 
can be circumvented in various ways, and one of them is construction of heterostructures, 
e.g. a semiconductor 'sandwiched' between a magnetic insulator and a superconductor. The 
similarity between the system Eq. (2.41) and the chiral p-wave superconductor is additionally 
supported by the presence of Majorana zero modes in vortices, when the condition (2.44) is 


satisfied. 


Nodal topological superconductivity 


The superconducting order parameter in many unconventional superconductors vanishes in 
special points (or sometimes lines) called ‘nodes’. As gapless systems, they are not classified by 


the ten symmetry classes mentioned in Sec. 2.2, however, they still host gapless, topologically 
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protected boundary modes. Additionally, the topological invariants that characterize nodal 
topological superconductors are momentum dependent. There are many systems predicted 
to realize this exotic phase, e.g. noncentrosymmetric superconductors [53] and Weyl super- 
conductors [54], however as is usual with intrinsic unconventional phases, it is difficult to 
unambigously prove that they exist in nature. It is therefore desired to construct a nodal 
topological phase from simpler ingredients, just like in the case of chiral superconductivity 
designed by the mixture of singlet pairing, magnetic field and Rashba interaction. One such 
example is a nodal topological state of Ising superconductors. A class of materials called 
transition metal dichalcogenides (TMD) exhibits the unconventional Ising pairing. Viewed 
from the out-of-plane direction, TMDs show a honeycomb lattice with broken sublattice 
symmetry. This y > —y mirror symmetry breaking results in the so called Ising spin-orbit 
coupling. In contrast to Rashba coupling, it pins the electron spins to out-of-plane directions, 
and because of that Ising superconductors are remarkably tolerant of in-plane magnetic 
fields, with critical fields exceeding the Pauli limit by up to 6 times [55]. One can therefore 
expect that the combination of s-wave superconductivity (present in TMDs), Zeeman field 
and spin-orbit interactions will lead to the emergence of a topologically non-trivial phase. 
Indeed this is the case, as shown by He et al. [56]. The Ising coupling vanishes along the 
high symmetry lines Il' — M in the hexagonal Brillouin zone. It is at the intersection of the 
normal Fermi surface and this line, that the superconducting gap vanishes upon application 
of the external, in-plane Zeeman field, leading to six pairs of point nodes with opposite 
chirality. Even though the time reversal symmetry is broken by the external field, an effective 
symmetry — a combination of time reversal and mirror symmetries is present, which together 
with the particle-hole symmetry ensure that there is also chiral symmetry. The system then 
falls into the BDI class (one of the ten symmetry classes mentioned before [39]) for any fixed 


k,, and can be characterized by the winding number 


W(k,) = f dk, Tr [CHp 3k, Hy] , (2.45) 


x 271 


with C the chiral symmetry operator. Whenever the system is in the state with W(k,) # 0, 


flat Majorana bands connect the point nodes. Since they are flat, they do not have any 
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velocity and are localized on the edges on the system, which are not parallel to the zigzag 
termination. When projected onto this direction, the nodal points cancel each other. The 
proposal for detection of such localized Majorana states is therefore scanning tunneling 
spectroscopy. Below the critical magnetic field we expect a usual s-wave gap when probing 
the edge, whereas above the critical field (in the topologically non trivial state), there would 


be a V-shaped gap (characteristic for nodal states) with a pronounced zero-bias maximum. 


26 


Chapter 3 


Methods 


3.1 Bogoliubov-de Gennes equations 


The initial applications of the Bogoliubov-de Gennes approach (BdG) consisted of analysis of 
quasiparticle exctitations around vortices. Significant increase in computing power allows 
for a treatment of more complicated structures and situations. The main appeal of the 
BdG method is the relative simplicity and attainability. Despite generically being a set of 
generalized Schródinger equations, BdG approach is widely used to tackle a rich variety of 
phenomena, such as the exotic Majorana quasiparticles in a chain of magnetic adatoms [57], 
disorder in d-wave superconductors [58], or unconventional superconductivity induced by 
proximity of d,, superconductor to a topological insulator [59], to name just a few. As 
mentioned above, the advances in local probing techniques, like the STM, revived the desire 
to theoretically access single-particle density of states in atomic scale, and for that the 
BdG approach is definitely apt. Below we will succinctly derive the Bogoliubov-de Gennes 
equations in their tight-binding version. This approach is generic, although we will essentially 
follow the program of Ref. [60], aside from slight differences in notation. We start with a 


Hamiltonian: 


= t cac! T T * 
H = ` hijoCiglio + `> Mii CigCjo! + V P» (Ach, + A; m) (3.1) 
ij) 0,0 i 


(ij)n,c (ij , 
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which is general enough to take into account various possibilities. The single particle term 
hijo contains spin-conserving hopping up to n-th nearest neighbors, but can also be equipped 
with different on-site terms, like the chemical potential (whose local change would serve as 
a non-magnetic impurity), sublattice potentials in appropriate systems and magnetic fields. 
Present is also a possible spin-orbit interaction, expressed as a hopping with a change in spin 
direction proportional to AU . The last term expresses the on-site pairing with the pairing 
potential V < 0 and order parameter A; = (cici). The structure of the Hamiltonian obliges 


us to mind both spin and particle-hole spaces, thus computing four commutators: 
cit, H] = » hijsch + 2, Aci + VAE, 

TX , 

[cm H| = — » hjach. — 2, AT cj — VA} Cj 


a (3.2) 
ils H| = 2, hag Cay + 2, AD Cj, — V Auch, 


eis] = 7 had SA + V Ate, 


instructs us that the Bogoliubov-Valatin transformation is a legitimate approach. We 
therefore express the electronic operators as linear combinations of electron-like and hole-like 


quasiparticles 


f 1 


Cia = X (Uon UU). Cis = D (UE — mU); (3.3) 


n n 


where the prime denotes summation over positive energy states labeled by n, and e = +1. 
We impose the effective Hamiltonian to be diagonal H, = >> Ey + Ens, and it is easy 
to check that . 

Lyn? , Hy] = (7) Ei. (3.4) 


We now substitute Eq. (3.3) into the commutators Eq. (3.2) and by comparing the coefficients 


of terms with 4,, and qÅ from Eq. (3.4), we obtain the BdG equations, written concisely in a 
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matrix form: 


»» io; = Enĝi, (3.5) 
j 
where 
hep a dO A 
, XL hny Ai 0 
AN aue " (3.6) 
0 A; hin ^ij 
and 
Ut 
^ Uig 
ĝi = . (3:7) 
Vit 
Vit 


In order to access the effects which various elements of the Hamiltonian have on the spatial 
structure of the order parameter, we need to solve the set of BdG equations with a self- 
consistency condition and diagonalize the Hamiltonian matrix Eq. (3.6), first with a random 


or guessed distribution of A;, and compute it in each iteration using 


1 


V n, n* n, n* E, 
A; = z 2 (us vi T Uj Vip )tanh (7) ; (3.8) 


n 


until a desired accuracy is achieved. In most cases we are interested in the local density 

of states, which is proportional to the imaginary part of the retarded Green’s function of 
1 

the system G; (E) via the formula pio (E) = ——Im(G;,(E)). Using the Bogoliubov-Valatin 
T 


transformation Eq. (3.3) in a familiar expression Gj, = ((cio|cl,}}) we readily arrive at 
pia (E) = \ting|?6(E — En) + Wino 6(E + En). (3.9) 


The Dirac delta is usually artificially broadened using Lorentzian or Gaussian distributions. 
We can see from the above reasoning, that after establishing the physical implications of a 
model we wish to study and thus constructing a Hamiltonian on a discretized lattice, our 


aim is to diagonalize the appropriate matrix and calculate the observables of interest using 
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the eigenvectors and eigenenergies. This task can be computationally demanding, as the size 
of matrix Eq. (3.6) is 4Nx4N with N the number of lattice sites. Simulating large lattices 
is therefore problematic and one usually resorts to using exaggerated values of A, since for 


smaller systems, the average level spacing can be of the order of the superconducting gap. 


3.2 Bond currents 


Acquisition of the eigenvectors and eigenvalues of Eq. (3.6) can be additionally used for 
calculation of the bond current, i.e. the average flow of charge between sites of the lattice. 
Their relation with Shiba states was first established in Ref. [61], where the authors have 
shown that through the magnetoelectric effect, the supercurrents around the impurity are 
carried by the YSR states. To obtain the expression for the current we start from the 


continuity equation 


ay DE Jig = 0, (3.10) 


J 


where p; = >> cl Ci is the density operator. We now use the Heisenberg equation 


oO 


th mt 


= [p H], wy 


where H is the appropriate Hamiltonian of the studied system. As we are interested in 
obtaining a vector field which visualizes the flow of charge between sites, we drop any term 
that induces on-site charge fluctuations. This process leaves only the hopping terms of the 
Hamiltonian contributing to the final expression. After the evaluation of commutators and 
insertion of the Bogoliubov-Valatin transformation Eq. (3.3), we obtain for the current from 
site 2 


J dXX uif (En) + ooun f(E) -ed. —— G2) 


j,0,0' 

Here we work in the natural units e = A = 1. In the above expression hee will depend on 
the type of the hopping. It can be a usual, spin conserving hopping between the neighbors of 
arbitrary rank, and/or various strains of spin-orbit interactions. The magnetoelectric effect 


in the context of magnetic structures on a superconducting surface is neatly visualized in 
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the following argument. The effective Hamiltonian of the surface of a three dimensional 
material will contain a Rashba term, stemming from the lack of inversion symmetry, i.e. 
Hy = k?/2m +2X(o x k) - 2. Therefore, an additional, spin-dependent contribution to the 
current will emerge, as can be seen by evaluating Oy Hj = k/m + AZ x ø. Therefore the 
question of such effect being present because of a single impurity is perfectly valid. One 
additional remark is due in the context of bound states around magnetic impurities. As 
mentioned before, the Shiba states carry the supercurrents, as in the language of the T' matrix 
approach, only the poles associated with them give rise to non-zero current. Furthermore the 
direction of the flow is reversed after the quantum phase transtition described in the previous 
Chapter. This statement is modified, and the situation more complicated when the structure 


of nee contains non-trivial terms (cf. description of paper III in Chapter 4). 


3.3 Majorana polarization 


Since the first report of the elusive Majorana quasiparticles being detected in a mesoscopic 
transport experiment |62], there has been a good amount of debate, whether the zero-bias 
conductance peaks are really Majoranas. Both the Shiba states and Andreev bound states 
can without a doubt exist at zero energy in special cases, while the Kondo resonance is 
guaranteed to be centered around the Fermi level. Theoretical studies are naturally helpful 
in such a situation and provide different interpretations and ways to test the validity of 
experimental reports. An example is the proposal of measuring the profile of the supercurrent 
in a junction consisting of two superconductors connected by a nanowire [63]. As one can 
never obtain a true zero energy Majorana bound sate when using numerical methods, and a 
state which is arbitrarily close to zero is not an indisputable proof of a topologically non-trivial 
state, it is useful to possess various computational tools at one’s disposal. When exploring 
systems harboring the exotic Majorana quasiparticles, one can resort to a quantity called 
Majorana polarization, It was introduced in Ref. [64] as a complex version of local density 
of states, suitable for the description of Majorana quasiparticles, but was valid only for a 
subset of systems. A generalization of this approach was provided in Ref. [65]. One can 


obtain it from the particle-hole operator, of which the Majorana bound states are eigenstates 
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with eigenvalue |1|. In the Nambu basis V; = (cis, ci; B» —cl,)? the particle-hole operator 
takes on the form C = e?g, T, K , with K as the complex conjugation operator. If we are 
interested in the local Majorana polarization vector, we evaluate, using a wavefunction 


Pi = (Uit, popa vi), the following: 
(IC) = -2) ousi (3.13) 


with C; = Cfi, where f; is the projection onto site ?. It is useful to inspect the real space map 
of Eq. (3.13), presented as a vector quantity, with the real and imaginary parts as r and y 
components. This way one can inspect if the state is a ’true’ Majorana, for to be considered 
as such, in a region R where we expect a Majorana state to be localized, the vectors must be 


'ferromagnetically' aligned. One can also verify this condition by evaluating whether 


x wie] 
iER Ee (3.14) 


» iV) ^^ 


iER 


when the region R is appropriately chosen. This method was succesfully employed to examine 
many different types of systems, as one dimensional wires, two dimensional strips or junctions 
between a normal metal and a superconductor [65,66]. It can as well be used to study the 
topological phase diagram of recently fabricated, highly controllable topological Josephson 
junctions (cf. the description of paper VI in Chapter 4). 
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Chapter 4 


Original contributions 


We will now present a summary of results obtained in the papers constituting the thesis and 
state the author's input to them. The first paper emerged from a wish to check if there exist 
mechanisms, which may influence the spatial character of the YSR wave function. Recent 
experimental results at the time revealed that the dimensionality certainly plays a crucial 
role in this context [67]. Seeing that in the lattice version, the spin-orbit interaction is a 
hopping term, we have hypothesized that it may influence the spatial extent of bound states. 
In the second publication we have summarized the efforts of our group and coworkers in 
the general context of bound states, whether coming from impurities, their collections or 
highly controllable systems — quantum dots. The Yu-Shiba-Rusinov states were presented 
in a tight-binding version for a generic situation of a magnetic impurity in a square lattice. 
Having learned that spin-orbit interactions affect the Shiba states in an unusual manner, we 
have turned our attention to different flavors of the coupling between orbital and spin degrees 
of freedom. It turned out that the interactions capable of inducing the exciting topologically 
non-trivial phases were an interesting example. We have therefore studied the effect of the 
intrinsic spin-orbit coupling of honeycomb lattices, introduced as a necessary ingredient for 
the presence of the quantum spin Hall efffect in graphene [9] on the bound states induced by 
single magnetic impurity. The next paper was the result of reflecting on recent experimental 
breakthroughs — discovery of two dimensional ferromagnetism |68] and Ising superconductivity 
in transition metal dichalcogenides [55,69]. We have envisioned a system, already known 


to harbor a rich topological phase diagram — a Shiba lattice [70,71], but with important 
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differences. This time, it was the combination of Ising spin-orbit coupling (bearing similarities 
to the Kane-Mele term, albeit of different origin) and in-plane ferromagnetism, that yielded a 
topologically non-trivial, nodal superconducting phase. Being familiar with the Ising spin-orbit 
interaction proved useful in the project with the experimental group led by Peter Liljeroth 
from Aalto University. They have manufactured state of the art heterostructures using 
molecular beam epitaxy. Islands of ferromagnetic CrBr3 on superconducting TMD — NbSes 
revealed zero bias peaks on their edges when probed with scanning tunneling spectroscopy. We 
have established a theoretical model showing that this designer structure realizes topological 
superconductivity with 1D chiral Majorana modes on the edges of the islands. The most 
recent work, connected to a highly controllable topological Josephson junction focuses on the 
Majorana polarization introduced in Chapter 3. We have shown that its absolute value can 
be probed by polarized scanning spectroscopy and examined the influence of an electrostatic 
impurity on the spatial structure of Majorana bound states induced in the junction. Below 


we discuss the selected papers in more detail. 


I. Yu-Shiba-Rusinov states of impurities in a triangular lattice of NbSes with spin-orbit 
coupling, Phys. Rev. B 96, 184425 (2017) 
A. Ptok, S. Gtodzik, T. Domariski 


Inspired by experimental results [67] showing that when the superconducting substrate 
is (quasi) two dimensional, the Shiba state wave function will spread significantly 
further in the lattice, we have examined a magnetic impurity embedded in a triangular 
lattice superconductor. Through numerical calculations we have shown that a spin-orbit 
interaction with its effective field acting in the plane of the substrate, the spatial length 
of YSR states can be additionally increased. We have also determined that this peculiar 
type of spin-orbit coupling will slightly affect the value of the critical coupling between 


the impurity and the superconducting host, at which the quantum phase transition 
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II. 


occurs. We have examined the spatial patterns of the YSR states (proportional to STM 
results), and by computing the displaced moving average of the YSR LDOS resolved 
the particle-hole oscillations measured in the aforementioned experiment. Lastly, we 
have presented the intereference between the YSR signals coming from magnetic dimers, 
and the influence of their relative position on such patterns. 

Author's contribution: Analytical & numerical calculations. Preparation of the 


manuscript. 


Interplay between pairing and correlations in spin-polarized bound states, Beilstein J. 
Nanotechnol. 2018, 9, 1370-1380 
S. Glodzik, A. Kobiałka, A. Gorczyca-Goraj, A. Ptok, G. Górski, M. M. Maska, T. 


Domański 


This paper reviews phenomena associated with the presence of magnetic adsorbates 
in contact with superconductors. From single impurities, through a chain of magnetic 
atoms, to an interpaly between the Kondo effect and Majorana modes in a setup with 
a quantum dot. The first part focuses on Yu-Shiba-Rusinov (YSR) states induced by 
a single magnetic impurity in a square lattice superconductor and shows the increase 
of their spatial extent due to an in-plane spin-orbit field. This is especially visible 
in the particle-hole oscillations, which show an increase of the spectral weight of the 
Shiba states away from the impurity sites, by almost an order of magnitude. Other 
characteristic features of YSR states are also present, albeit their topography is naturally 
different than in the case of triangular lattice. 

Author's contribution: Analytical & numerical calculations presented in the para- 


graph ,Single magnetic impurity". 
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III. In-gap states of magnetic impurity in quantum spin Hall insulator proximitized to a 
superconductor, J. Phys.: Condens. Matter 32 (2020) 235501 
S. Glodzik, T. Domański 


In this paper we have studied a single magnetic impurity embedded in a honeycomb 
lattice with proximity induced superconducting order. Because of the symmetries in 
such a lattice, a spin conserving, next-nearest-neighbor hopping term, called intrinsic 
spin-orbit coupling, or the Kane-Mele term is allowed [9]. The presence of this coupling 
drives the system into the topologically non-trivial state — the quantum spin Hall phase, 
with a Z5 topological invariant [41]. We have established that the presence of the Kane- 
Mele term shifts the critical magnetic interaction Jc to higher values, and significantly 
reduces the spatial extent of Shiba-like states. T'hose are not strictly Yu-Shiba-Rusinov 
states, because the density of states in a honeycomb lattice vanishes at the Fermi energy, 
and the formula Eq. (2.27) would not yield a subgap state. Nevertheless there are 
bound states coming from the magnetic impurity, and we have observed that when the 
quantum spin Hall phase is not destroyed (e.g. by the chemical potential shift), there 
are two species of impurity bound states in the gap. One pair (Shiba-like) undergoes the 
quantum phase transition, while the other does not. There are usual manifestations of 
the change in the ground state: sign reversal of the order parameter at the impurity site, 
change in the bulk polarization and the reversal of the bond current. When inspected 
carefully, it turns out that only the current carried by the YSR-like states is reversed. 
Because of the competition of two flow directions, total current is drastically reduced 
after the phase transition. 

Author's contribution: Analytical & numerical calculations. Preparation of the 


manuscript. 
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IV. Engineering nodal topological phases in Ising superconductors by magnetic superstruc- 
tures, New J. Phys. 22 (2020) 013022 
S. Glodzik, T. Ojanen 


Motivated by the discovery of Ising superconductivity in TMDs [55,69] and the prediction 
of flat Majorana bands realized in this type of materials under in-plane magnetic 
fields |56], we have envisioned an island of classical spins pointing parallel to the plane 
on the surface of a superconducting TMD. Such a system is expected to enter a nodal 
(gapless) topological phase in which the flat Majorana bands would manifest as localized 
zero energy states. T'hose states are only present on the island edges parallel to the 
armchair termination of the lattice. This stems from the the structure of the point 
nodes (gap closings) in the momentum space. Their configuration additionally controls 
the degeneracy of the flat band, reflected in the value of the topological invariant — the 
integer winding number. 

Author's contribution: Analytical & numerical calculations. Preparation of the 


manuscript. 


V. Topological superconductivity in a designer ferromagnet-superconductor van der Waals 
heterostructure, arXiv:2002.02141 
S.Kezilebieke, M. Nurul Huda, V. Vaňo, M. Aapro, S. C. Ganguli, O. J. Silveira, 
S. Glodzik, A. S. Foster, T. Ojanen, P.Liljeroth 


Two dimensional topological superconductivity induced by magnetic impurities has 
been proposed theoretically [70,71] and observed experimentally |72, 73], altough there 
is some debate concerning the robustness and scalability of the experimental results. 
By designing a structure comprising two types of van der Waals materials — NbSes 
and CrBr3, which exhibit two competing electronic orders (superconductivity and 


ferromagnetism respectively), we show that the system is in the topological supercon- 
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VI. 


ducting state with pronounced zero-bias peaks on the edges of the islands. Those edge 
states, revealed by the STM data and theoretical calculations, are the celebrated chiral 
Majorana states, proposed as building blocks for topological quantum computation. 
The theoretical model correctly predicts that in the phase with the Chern number 
C — 3 the edge states are the most prominent excitation even at energies close to the 
topological gap. Their exponential localization at the island edges is also consistent, 
as well as the apparent (but not real) discontinuity of the edge state, stemming from 
interference effects due to irregularities of the boundary. 

Author's contribution: Establishing the theoretical model and performing the mo- 


mentum & real space numerical calculations. 


How to measure a Majorana: The Majorana polarization of a topological planar Joseph- 
son junction , arXiv:2004.01420 
S. Glodzik, N. Sedlmayr, T. Domański 


It has been first predicted theoretically [74,75] and then shown experimentally [76,77] 
that a topological Josephson junction can be driven into a topologically non-trivial 
state using the phase difference between the superconducting electrodes. The quasi- 
one-dimensional barrier between the two leads hosts Majorana bound states which are 
remarkably robust. We test this robustness by calculating the Majorana polarization for 
different widths of the barrier and show that this value can be probed using polarized 
Andreev spectroscopy. Additionally we examine the influence of a strong point-like 
electrostatic impurity and show that when placed close to the end of the normal 
region, it can localize one of the Majorana bound states, leaving the other one virtually 
unaffected. We also confirm the phase diagram using the polarization criterion. 


Author's contribution: Numerical calculations. Preparation of the manuscript. 
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Yu-Shiba-Rusinov states of impurities in a triangular lattice of NbSe; with spin-orbit coupling 


Andrzej Ptok,?:* Szczepan Glodzik,^ and Tadeusz Domafiski^ 
Vinstitute of Nuclear Physics, Polish Academy of Sciences, ul. E. Radzikowskiego 152, 31-342 Kraków, Poland 
? Institute of Physics, Marie Curie-Sklodowska University, pl. Marii Sklodowskiej-Curie 1, PL-20031 Lublin, Poland 
(Received 5 July 2017; revised manuscript received 31 October 2017; published 22 November 2017) 


We study the topography of the spin-polarized bound states of magnetic impurities embedded in a triangular 
lattice of a superconducting host. Such states have been observed experimentally in 2H-NbSe, crystal 
[G. C. Ménard et al., Nat. Phys. 11, 1013 (2015)], and they revealed oscillating particle-hole asymmetry extending 
to tens of nanometers. Using the Bogoliubov-de Gennes approach, we explore the Yu-Shiba-Rusinov states in 
the presence of spin-orbit interaction. We also study the bound states of double impurities for several relative 


positions in a triangular lattice. 


DOI: 10.1103/PhysRevB.96.184425 


I. INTRODUCTION 


Magnetic impurities are detrimental for Cooper pairs 
because they induce spin-polarized subgap states [1—3] and 
(when impurities are dense enough) partly suppress or 
completely fill in the energy gap of the superconducting 
sample. Such in-gap quasiparticles, dubbed Yu-Shiba-Rusinov 
(YSR) states [4—6], have been observed experimentally in 
various systems [7-14]. They always exist in pairs, appearing 
symmetrically with respect to the chemical potential (treated 
here as the "zero-energy" reference level). Their energies 
can be controlled either electrostatically or magnetically [14]. 
Another feature is their spin-polarization evidenced by the 
asymmetric conductance at opposite voltages [15,16]. 

The bound states formed at magnetic impurities in three- 
dimensional isotropic superconductors are usually character- 
ized by a relatively short spatial extent [17]. Contrary to that, 
in two-dimensional (2D) lattices, Ménard et al. [11] have 
reported different behavior, displaying a much longer extent 
of the YSR states with alternating (oscillating) particle-hole 
spectral weights. Furthermore, the bound states of impurities 
in superconducting 2H-NbSe; [18] with quasi-2D triangular 
lattice structure and strong spin-orbit coupling [19] have 
revealed the long-range coherent structures of a star-shape, 
whereas molecular dimers developed more complex spatial 
features [20]. 

Bulk crystals of 2H-NbSez are characterized by cen- 
trosymmetric (P63/mmc) structure, formed by the stacking 
of noncentrosymmetric layers [21,22] [Fig. 1(a)]. Every layer 
is arranged in a honeycomb structure, comprising Nb and Se 
sublattices. Local inversion symmetry breaking [23—25] gives 
rise to the out-of-plane spin polarization [21] in every layer. At 
Tcpw © 33 K, there appears charge-density order [22,26,27], 
and below T, ~ 7 K [21] the superconducting state sets in. 

The normal state Fermi surface, studied by angle-resolved 
photoemission spectroscopy (ARPES) [21,27—30], has re- 
vealed two pairs of Nb-derived pockets, which are trigonally 
warped around a central I’ point and at the corner of the 
(hexagonal) first Brillouin zone. Ab initio (DFT) calculations 
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indicated that the Fermi surface sheets originate predominantly 
from Nb 4d orbitals [21,31—33]. Asa consequence, a triangular 
lattice formed by Nb atoms plays an important role for the 
superconducting state of this compound [34] and implies 
further a unique starlike shape of the bound states. 

Contrary to bulk systems, the Fermi surface of the single 
monolayer 2H-NbSe; consists of only the pockets around the 
corner points of the Brillouin zone [21], whose size depends 
on the spin polarization [Fig. 1(b)]. The latter effect originates 
from the in-plane spin-orbit field [19,30]. Coupling of the spin 
to the valley distinguishes between nonequivalent parts of the 
Brillouin zone. Similar behavior has also been observed in 
other materials with hexagonal lattice structures [36—39]. 

Motivated by the recent experimental results of Ménard 
et al. [11], we study the YSR states of magnetic impurities 
embedded in a triangular lattice of the 2D superconducting 
host. The single monolayer of the 2H-NbSe» can be treated 
as a two-dimensional triangular lattice [34] with an in-plane 
spin-orbit field (in the supplemental material, we take into 
account also the out-of-plane spin-orbit component, which 
might be realized in other compounds [40,41 ]). 

In Sec. II we present the microscopic model and discuss 
some methodological details. Next, in Sec. I A, we analyze 
the YSR bound states of a single magnetic impurity in a 


FIG. 1. (a) Crystallographic structure of the centrosymmetric 
NbSe; compound and its primitive unit cell (black prism). The image 
was obtained using VESTA software [35]. (b) Schematic view of the 
Fermi surface in the NbSe; monolayer, adopted from Ref. [21]. Blue 
(red) corresponds to different spin [negative (positive)] polarizations 
for each Fermi sheet in the zone corners. 
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FIG. 2. Schematic illustration of a single magnetic impurity 
(yellow) in a two-dimensional triangular lattice. Each lattice site is 
surrounded by six nearest neighbors at positions dg. 


triangular lattice (Fig. 2), focusing on the role of spin-orbit 
coupling. In Sec. IIIB we study the bound states of double 
magnetic impurities (arranged in three different configura- 
tions) that might be relevant to the experimental data [20] 
revealing strong interference effects. Finally, in Sec. IV we 
summarize the main results. 


Il. MODEL AND METHOD 


Magnetic impurities embedded in a 2D superconducting 
host with spin-orbit coupling can be described by the following 
Hamiltonian: 


H = Ho + Himp + Him + Asoc. (1) 
The single-particle term 
Ho 2 -t Y de n 86s Q) 
li, jjo io 


describes the kinetic energy of electrons, where al (Gio) 
denotes the creation (annihilation) of an electron with spin 
c at the ith lattice site, £ is a hopping integral between the 
nearest neighbors, and jz is the chemical potential. Large spin 
S of the impurities can be treated classically [1,2], and in such 
a case the interaction potential can comprise the magnetic J 
and nonmagnetic K parts, 


Pimp = —I(C, 204 — êb ĉo) + KG, êo +6), C0). — 


We describe the superconducting state, imposing the on-site 
interaction 


Hint = U b3 0 ên (4) 
i 


with attractive potential U <0. Such effective pairing is 
assumed to be weak, therefore we can treat it within the 
standard mean-field decoupling 


SUNT. êi = xê et, + xféí63 — lxil 


+nnêl ên + ni Gen —nini, (5 
where x; = (6;,6;4) is the superconducting order parameter 
and nj; = (a, Ĉio) is the average number of spin o particles 
at the ith site. The Hartree term can be incorporated into 
the effective local spin-dependent chemical potential u > 
[lic = u — Unis. AS we shall see, impurities suppress the 


local order parameter x; whose magnitude and sign depend on 
the coupling strength J [42,43]. 
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The spin-orbit coupling (SOC) can be expressed by [44] 


Hsoc = —ix X` CNMCE x 67"). Ô ĉion, (6) 


ijooc' 


where the vector d; — (47 ,d;,0) corresponds to the nearest 
neighbors of the ith site (Fig. 2), and 6 = (0,,0,,0;) consists 
of the Pauli matrices. The unit vector 1? shows a direction of 
the spin-orbit field, which in general can be arbitrary, but we 
restrict our considerations to the in-plane ® = £ = (1,0,0) and 
out-of-plane 2 — (0,0,1) polarizations, so formally we have 


for in-plane field, 


x » 
dio, — do, 


j for out-of-plane field. 


(7) 


Let us notice that the out-of-plane component mixes 7 and 
| particles, whereas the in-plane field corresponds to addi- 
tional spin-conserving hopping with the spin- and direction- 
dependent amplitude. 


Bogoliubov-de Gennes technique 


Magnetic impurities break the translational invariance of 
the system, therefore the local pairing amplitude x; and 
occupancy Mio have to be determined for each lattice site 
individually [45]. One can diagonalize the Hamiltonian (1) 
via the following unitary transformation: 


ĉio = Y (ino Pa — 0v). (8) 


where ?,, and i are quasiparticle fermionic operators, with the 
eigenvectors Ujng and Vino. This leads to the Bogoliubov-de 
Gennes (BdG) equations [46] 


Uint ij U jnt 
efom] Y Dj -H& 0 Sj || any 
n = 
Uin} 7 ap 0 Hij, Dij Ujin} 
; N j 
Vint 0 5i; Di; -Hý Ujnt 
(9) 
containing the single-particle term Ajj, = —19(i,j; — 
[is + (K — o J)8io]ó;; + Sj? and the spin-orbit coupling 
ien gue = —iAY, (dı x 67" y. db ój;4q,. Here, S7? and 


SP? correspond to the in-plane and out-of-plane spin-orbit 
field, respectively, which satisfy Spe = (S77 )", and D;; = 
Ux;6;; describes the on-site pairing. The superconducting 
order parameter x; and occupancy n;, can be computed 
self-consistently from the BdG equations (9), 


Xi = 3 s v f) — Ming Uy ACE O) 


nio = ) linc f(&) + [vina PCE), — (10 


n 


where f (c) = 1/[1 + exp(w/kgT)] is the Fermi-Dirac distri- 
bution. In particular, the spin-resolved local density of states 
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(LDOS) is given by [47] 


pio (e) = Y Misc 0 — En) + [ine ?8(@ + En). (12) 


n 


For its numerical determination, we have replaced the Dirac 
ô function by the Lorentzian 5(w) = Z/[r (e? + ¢7)] with a 
small broadening ¢ = 0.025t. 


II. NUMERICAL RESULTS AND DISCUSSION 


We now present the BdG results obtained for the single 
impurity embedded in a triangular lattice (Sec. III A) and for 
several configurations of two magnetic impurities (Sec. III B). 
Numerical computations have been done at zero temperature 
for the finite cluster Na x N, = 41 x 41, assuming U/t = 
—3, w/t = 0, K/t = 0, and determining the bound states for 
varying J. In this work, we focus on the effect of an in-plane 
spin-orbit field, and additional results for the out-of-plane SOC 
are shown in the Supplemental Material (SM) [48]. 


A. Single magnetic impurity 


Let us start by discussing the results obtained in the absence 
of spin-orbit coupling. A typical quasiparticle spectrum is 
displayed in Fig. 3, where we can recognize the gapped 
region |w| € A of the superconducting host (for our set of 
model parameters A ~ 0.651) and one pair of the bound states, 
appearing symmetrically around the chemical potential. The 
energies +E, of these states and spectral weights depend on 
the coupling J. In particular, at some critical J, (indicated by 
black arrows) they eventually cross each other. This crossing 
is a hallmark of the quantum phase transition (QPT) [49] 
in which the ground state undergoes a qualitative evolution 
[15]. When magnetic coupling overcomes the pairing energy 
(for |J| > Je), the particle and hole states become degenerate, 
and the ground state changes from a BCS singlet to a spinful 
configuration [15,50—52]. 


"5.0 -2.5 0.0 2.5 5.0 
J/t 


FIG. 3. Evolution of the low-energy spectrum with respect to the 
magnetic coupling J. The solid black line shows the magnitude of 
the pairing gap in regions far away from the impurity, black arrows 
point to the quantum phase transition (i.e., crossing of the subgap 
YSR states), and thin-dashed lines display the YSR bound states 
calculated from Eq. (13). Results were obtained without the SOC. 
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FIG. 4. Influence of the in-plane SOC on the critical value Je. 
Blue and red correspond to the discontinuous change of U xo at the 
impurity site, and a white line marks the QPT. Black arrows indicate 
two values of 4, for which the profiles are shown in Fig. 5. 


Our BdG data can be confronted with the analytical results 
of the thermodynamic limit Na x N, — oo [50]: 
fecu = 13 
ysr ar IL ni (13) 
where o = zpoJ is the dimensionless impurity coupling 
parameter, po is the normal state DOS at the Fermi level, and 
A is the superconducting gap. These quasiparticle energies 
(13) are displayed in Fig. 3 by a thin-dashed line. In the weak- 
coupling limit |J| < Je, the formula (13) matches well with 
our numerical BdG results. Some differences appear above 
the QPT (for |J| 2 Je), where the local pairing parameter 
at magnetic impurity is substantially reduced, affecting also 
the pairing gap of its neighboring sites. With an increasing 
coupling A, the QPT is shifted to higher values (Fig. 4). The 
critical Je corresponds to a value of J at which the YSR states 
cross each other. Variation of the critical Je is caused by the 
influence of the SOC merely on the normal state DOS (jo). 
Such a quantum phase transition is manifested by a sign 
change ofthe order parameter xo at the impurity site [Fig. 5(a)], 


(a) o3 


02r 


0.1 


Xo 


0.0 25 5.0 


J/t 


FIG. 5. The order parameter xy = (69,694) obtained at the impu- 
rity site i = 0 (a) for the weak (red line) and strong (blue dotted line) 
spin-orbit couplings, with A/t = 0.1 and 1.0, respectively. Magnetic 
polarization of the YSR states ~9;(@) — po, (cw) (b), obtained for 
A/t — 0.1. 
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FIG. 6. Spatial patterns of the “negative” and “positive” YSR 
states fos + Bil Results are obtained for J/t = +1.2 (top panel) 
and J/t = +2.5 (bottom panel), assuming the in-plane spin-orbit 
coupling A/t = 0.1. 


and discontinuity of its absolute value is a signature of the 
first-order phase transition [53-55]. Let us emphasize that the 
QPT is associated also with a reversal of the YSR polarization 
[Fig. 5(b)] and furthermore the total polarization of the system 
P= 1 Y^; (it — ni,) abruptly changes at J = +J. from zero 
to +1/2 [56]. Similar behavior can be observed for multiple 
impurities [57]. 

In the weak-coupling limit (i.e, for A « t), we can 
hardly notice any meaningful influence of SOC on the bulk 
superconductivity and the YSR states (see Fig. 1 in the SM 
[48]). A similar conclusion has been previously reported from 
the T-matrix treatment of magnetic impurities for ID and 2D 
square lattices by Kaladzhyan et al. [52]. Our calculations 
have been done for 4/t = 0.1 which could be realistic for 
the NbSe; compound. Obviously for much stronger values 
of the spin-orbit coupling, both the superconducting state and 
the bound YSR states depend on the magnitude of A and the 
direction of the magnetic moment [58]. 

Let us now explore the spatial extent of the YSR states. 
This can be achieved within the BdG approach by integrating 
the spectral weights 


02 
p= f podo (14) 
on 
in the interval w € (@;,w2) capturing every in-gap quasipar- 
ticle below or above the Fermi level [59]. Our numerical 
calculations have been done for the single impurity in the 
weak J = —1.2t > — J. and strong magnetic coupling limits 
J = —2.5t < — Je, respectively. The results shown in Fig. 6 
(notice different scales for each of these panels) reveal the 
characteristic six-leg star shape, whose extent spreads on 
several sites around the magnetic impurity. Spectral weights at 


PHYSICAL REVIEW B 96, 184425 (2017) 


J/t=+1.2 
| — NG | =g" 
d ler + enll ol meee less + Pali 
: i Jes 
E 0 
35 
T 
JE = 
0 


FIG. 7. The same as in Fig. 6, but for the stronger spin-orbit 
interaction A/t = 1. 


the positive and negative energies are quite different, leading 
to a finite spin polarization of the YSR states (displayed in the 
bottom panel in Fig. 5). 

Upon varying J, we observe that the star-shape (charac- 
terizing the Cg symmetry of a triangular lattice) is rather 
robust. Such YSR state patterns could be probed by scanning 
tunneling microscopy, which nowadays has an atomic-scale 
resolution [8,12,60—62]. Let us emphasize that this six-leg 
star shape originates from a triangular lattice geometry and 
from the particular topology of the Fermi surface [63], in 
agreement with the experimental observations [11]. We have 
checked that the YSR states are only quantitatively (by a few 
percent) affected by the weak (4./t = 0.1) in-plane spin-orbit 
coupling. For stronger SOC, the results are presented in Fig. 7. 
Comparison with Fig. 6 shows that for higher A the YSR 
states gradually increase their extent, and the starlike shape 
seems to be weakly deformed with some elongation parallel 
to the x axis. Figure 7(a) presents an especially large range 
of bound states, extending well beyond ten lattice constants 
from the impurity site. As stated in the previous section, the 
in-plane spin-orbit field leads to the presence of the 57^ term 
in the single-particle part of the BdG equations, which directly 
affects the hopping amplitude. To observe the influence of this 
term on the spectral function, A has to be of the order of t. We 
suspect that the large spatial extent of the YSR states reported 
in [11] was a consequence of the reduced dimensionality 
and/or the structure of the atomic lattice, and was not affected 
much by the in-plane spin-orbit field of NbSe». In general, 
however, materials with the stronger SOC couplings could 
reveal some increase in the spatial extent of subgap bound 
states. 

For some quantitative analysis of the spatial profiles of 
YSR states, we define the displaced moving average (DMA) 
p^ (r) interpreted as an averaged spectral weight contained in 
a ring of radius r and its half-width ór. It depends only on 
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FIG. 8. Profiles of the hole- (blue line) and electron-like (red line) 
displaced moving average (DMA) for the YSR bound states p*(r) as 
a function of distance r from the impurity (with dr = 0.5). The left 
and right panels correspond to |J| < Je and |J| > Je, respectively. 
Results are obtained for the in-plane spin-orbit field A/t = 0.1. The 
dashed gray line corresponds to exp(—r), which is a guide to eye. The 
blue (red) background indicates the dominant hole (particle) type of 
the YSR state at a given r. 


a radial distance r from the magnetic impurity ro, averaging 
the angle-dependent fluctuations. Our results are presented 
in Fig. 8. They clearly show that functions p~(r) of the 
YSR states are characterized by particle and hole oscillations 
that are opposite in phase (see the blue and red lines). Such 
particle-hole oscillations decay exponentially with distance 
(notice a logarithmic scale in Fig. 8), in agreement with 
previous studies [11,56,64]. In the 2D continuum version of 
this model, the wave functions of the YSR states have been 
expressed analytically [11]: 


w*(r) x ET. sin (ter = 1 x s) 


7 r 
x exp |- sin(ó* — DE| (15) 


where kpr is the Fermi wave vector, r is the distance from the 
impurity, whereas ¢ is the superconducting coherence length. 
Both functions oscillate with krr, but with different scattering 
phase shifts 5+ = Kpy + J/po. At short distances, the YSR 
wave functions are governed by sin(krr)/./krr, whereas for 
larger r the exponential envelope function suppresses particle- 
hole oscillations (see the dotted line in Fig. 8). Dominant 
(particle or hole) contributions to the YSR bound states are 
displayed by an alternating color of the background in Fig. 8. 
The period of such oscillations is approximately equal to nearly 
two lattice constants. The out-of-plane spin-orbit field leads to 
a similar behavior (Fig. 5 in [48]). 

The quantum phase transition (at Je) has consequences on 
a reversal of the magnetization induced near the impurity (see 
Fig. 9). For |J| < Je, the impurity is weakly screened, whereas 
for stronger couplings, |J| > Je, the impurity polarizes its 
neighborhood in the direction of its own magnetic moment. In 
both cases, this short-range magnetization does not coincide 
with the six-leg star shape of the bound states. Differences 
between the YSR wave functions and various components of 
magnetization have been previously discussed for a 2D square 
lattice by Kaladzhyan et al. [52]. 
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FIG. 9. Magnetization along the z axis induced near the magnetic 
impurity for |J| < Je (a) and |J| > Je (b). 


B. YSR of double impurities 


The BdG technique has the advantage that it can be easily 
applied to study the bound states of more numerous impurities, 
distributed at arbitrary positions in a crystallographic lattice. 
In this section, we consider the case of double magnetic 
impurities arranged in three different configurations displayed 
in Fig. 10(a). Our study of the YSR states is inspired 
by the results of Ref. [20] for ferromagnetic dimers. Such 
BdG calculations can be applied to more complex molecules 
[51,65-67] and/or multi-impurity structures [68—70]. It is 
well known [1] that multiple impurities can develop several 
quantum phase transitions with some characteristic features. 
They have been previously studied for 2D lattices, treating 
the spins classically [57,71,72] and taking into account the 
strong correlation effects within the Anderson-type scenario 
[73]. Here we explore the YSR states of two classical magnetic 
impurities embedded in a triangular lattice, assuming the weak 
in-plane spin orbit interaction A/t = 0.1. 

Figure 10 presents the subgap spectrum obtained for 
different configurations of the double impurities. We notice 
that coupling between the impurities induces the double-peak 
structure of YSR states, both at negative and positive energies 
[panel (b)]. Figure 11 displays spatial distributions of the 
YSR states for each configuration of the double magnetic 
impurities for the weak (left column) and strong (right column) 
couplings J. Although the Cg rotational symmetry is broken, 
one can clearly see the mirror symmetry with respect to the axis 


(a) (b) 
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FIG. 10. Schematic illustration of two magnetic impurities ar- 
ranged in three different configurations (a) and the subgap spectrum 
(b) for the nearest neighbors (NN), next nearest neighbors (NNN), and 
third nearest neighbors (3NN), as shown by red, blue, and green lines, 
respectively. We assumed the in-plane spin-orbit coupling A/t = 0.1. 


184425-5 


PTOK, GLODZIK, AND DOMANSKI 


J/t = 1.2 z J/t = 2.5 
— NN | 
140, 60 
Mi i oL 
[| 
[| 
‘po Qe 
d "u EAE JE 
a pay 
[| 
40L- 4 Bob 4 
[| 
1o D 10 O 4 40 0 16 0 
N 


FIG. 11. Spatial pattern ofthe YSR states for the double magnetic 
impurities in different (NN, NNN, 3NN) configurations, as indicated. 
The left column panels refer to the weak coupling J/t = —1.2 and 
the right column panels to the strong coupling J/t = —2.5 limits, 
respectively. We imposed the in-plane spin-orbit coupling à/t = 0.1 
that could be relevant to NbSe;. 


connecting these double impurities and the axis perpendicular 
to it. Novel spatial patterns of the YSR states are due to 
the constructive/destructive quantum interference between the 
overlapping subgap states. Obviously, the most significant 
quantum interference occurs for the quantum impurities either 
at the nearest-neighbor (NN) or next-nearest-neighbor (NNN) 
configurations, with a clear bonding-antibonding splitting of 
the YSR quasiparticle energies. For more distant arrangements 
of the double impurities (for instance, 3NN), their spatial 
patterns gradually evolve back to the starlike shape. A more in- 
depth comparison ofthe results obtained with and without SOC 
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is presented in [48] Figs. 6 and 7 of the SM. We hope that our 
theoretical predictions can be empirically verified, using the 
combined AFM (capable of manipulating the impurities) and 
STM (suitable for probing the subgap spectrum) techniques. 


IV. CONCLUSIONS 


In summary, we have investigated the energies and spatial 
extent of the Yu-Shiba-Rusinov (YSR) states induced by the 
classical magnetic impurities embedded in a two-dimensional 
triangular lattice of the 2H-NbSe; superconducting host. To 
study this particular crystallographic geometry in the presence 
of local inhomogeneities (in a form of the single or double 
impurities) and the spin-orbit coupling (SOC), we have 
adopted the Bogoliubov-de Gennes formalism. 

In agreement with experimental observations [11], we have 
found that the YSR states acquire sixfold rotational symmetry 
(starlike patterns) whose spectroscopic signatures extend onto 
about a dozen of the intersite distances. Furthermore, their 
intertwining (z -shifted in phase) particle-hole oscillations are 
clearly visible. The weak spin-orbit coupling (which should be 
relevant to the 2H-NbSez compound) has a rather negligible 
influence on the energies of YSR states, but for relatively 
stronger SOC their spatial extent eventually increases (beyond 
10 lattice constants in some cases). Analysis of the SOC for 
the single impurity indicates that the extended range of YSR 
states reported in [11] stems from the dimensionality and/or 
the structure of an atomic lattice, rather than from the in-plane 
spin-orbit field of such materials. 

We have also studied the subgap quasiparticle spectrum of 
double magnetic impurities in three different configurations, 
revealing either the constructive or destructive quantum 
interference that breaks Cg symmetry of the YSR wave 
functions. Deviation from the starlike shape (typical for single 
impurities) depends on the relative distance between such 
magnetic impurities. When they are close to each other, the 
YSR states develop a double-peak structure (characteristic for 
the bonding and antibonding states) whose spatial patterns no 
longer resemble the star shape. With an increasing distance 
between the impurities, such bonding-antibonding splitting 
gradually disappears, and the spatial starlike shape of YSR 
states is gradually restored. 
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We investigate single and multiple defects embedded in a superconducting host, studying the interplay between the proximity-in- 


duced pairing and interactions. We explore the influence of the spin-orbit coupling on energies, polarization and spatial patterns of 


the bound (Yu-Shiba-Rusinov) states of magnetic impurities in a two-dimensional square lattice. We also address the peculiar 


bound states in the proximitized Rashba chain, resembling the Majorana quasiparticles, focusing on their magnetic polarization that 


has been recently reported by S. Jeon et al. (Science 2017, 358, 772). Finally, we study leakage of these polarized Majorana quasi- 


particles into side-attached nanoscopic regions and confront them with the subgap Kondo effect near to the singlet-doublet phase 


transition. 


Introduction 

Magnetism is usually detrimental to superconductivity because 
it breaks the Cooper pairs (at the critical field strength H¿2). 
There are, however, a few exceptions in which these phenome- 
na coexist, e.g., in iron pnictides [1], CeCoIns [2]. Also, some- 
times magnetic fields induce superconductivity [3]. Plenty of 
other interesting examples can be found in nanoscopic systems, 
where magnetic impurities (dots) exhibit a more subtle relation- 
ship with the electron pairing driven by the proximity effect 
[4.5]. Cooper pairs easily penetrate the nanoscopic impurities, 
inducing the bound (Yu-Shiba-Rusinov) states that manifest 


the local pairing in coexistence with magnetic polarization. 
Such bound states have been observed in various systems 
[6-14]. In-gap states (appearing in pairs symmetrically around 
the Fermi level) can be nowadays controlled electrostatically or 
magnetically [12] whereas their topography, spatial extent and 
polarization can be precisely inspected by the state-of-art 


tunneling measurements [15,16]. 


It has been reported that adatoms deposited on a two- 
dimensional (2D) superconducting surface develop 
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Yu-Shiba-Rusinov (YSR) states, extending to a dozen of inter- 
site distances and they reveal particle-hole oscillations [11]. 
Bound states of these magnetic impurities in superconducting 
NbSe; are characterized by the star shape [17] typical for the 
rotational symmetry of its triangular lattice. More complex 
objects, such as dimers, reveal other spatial features, showing 
the bonding and antibonding states [18]. In a somewhat differ- 
ent context it has been pointed out [19] that exchange coupling 
between numerous quantum defects involving their intrinsic 
spins can couple them ferromagnetically. This can be used (e.g., 
in metallic carbon nanotubes) for a robust transmission of mag- 


netic information over large distances. 


In all cases the bound YSR states are also sensitive to interac- 
tions. One of them is the spin—orbit coupling (usually mean- 
ingful at boundaries, e.g., surfaces) [20-22]. Such interaction in 
one-dimensional magnetic nanowires can induce the topologi- 
cally nontrivial superconducting phase, in which the YSR states 
undergo mutation to Majorana (zero-energy) quasiparticles. 
Coulomb repulsion between the opposite spin electrons can 
bring additional important effects. In the proximitized quantum 
dots it can lead to a parity change (quantum phase transition) 
with further influence on the subgap Kondo effect (driven by 
effective spin-exchange coupling with mobile electrons). 
Furthermore, such spin exchange can be amplified by the in- 
duced electron pairing, and can have constructive influence on 
the Kondo effect [23,24]. 


We study here the polarized bound states, taking into account 
the spin-orbit and/or Coulomb interactions. In particular, we 
consider: (1) a single magnetic impurity in a 2D square lattice of 
a superconducting host, (ii) a nanoscopic chain of magnetic 
impurities on the classical superconductor (i.e., proximitized 
Rashba nanowire) in its topologically trivial/nontrivial super- 
conducting phase, and (iii) a strongly correlated quantum dot 
side-attached to the Rashba chain, where the Kondo and the 
leaking Majorana quasiparticle can be confronted with each 
other. These magnetically polarized YSR and Majorana quasi- 
particles as well as the subgap Kondo effect can be experimen- 
tally verified using tunneling heterostructures with ferromag- 
netic lead (STM tip). 


Results and Discussion 

Single magnetic impurity 

Let us start by considering a single magnetic impurity on the 
surface of an s-wave superconductor in presence of spin—orbit 
interactions. This situation can be modeled by the Anderson- 
type Hamiltonian 
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We describe the superconducting substrate by 


Hsc =—t » Gey + UD hyn, — 03 fias Q) 
i 


(ij) io 


where êl (ĉis) denotes creation (annihilation) of an electron 
with spin o at the i-th site, t is a hopping integral between the 
nearest neighbors, u is the chemical potential, and 7j; = ae 
is the number operator. For simplicity, we assume a weak 
attractive potential U < 0 between itinerant electrons and treat it 
within the mean-field decoupling 

ése ên = êh et exiesg en - Du P 


enge T ngêt ên TAM; 


where x; = (6,64) is the local superconducting order parame- 
ter and nj, = (nj). The Hartree term can be incorporated 
into the local (spin-dependent) chemical potential 
U> fg = — Ung. The second term in Equation 1 refers to 
the local impurity 


Pimp =-J (Ey it io +K (Argon iL). 0 


which affects the order parameter y; near the impurity site i = 0, 
inducing the YSR states [25,26]. In this work we focus on the 
magnetic term J [4,27], disregarding the potential scattering K. 


The spin—orbit coupling (SOC) can be expressed by 


where the vector d; = (dj, d7,0) refers to positions of the 
nearest neighbors of the i-th site, and © = (Ox, Oy, 0;) stands for 
the Pauli matrices. The unit vector w shows the direction of the 
spin-orbit field, which can be arbitrary. Here we restrict our 
considerations to the in-plane wzií- (1, 0, 0) polarization, 
which will be important for nontrivial superconductivity in 
nanowires discussed in the subsection 'Magnetically polarized 
Majorana quasiparticles’. The other (out-of-plane) component 
could eventually mix 1 and | spins [22]. 


Impurities break the translational invariance, therefore the 
pairing amplitude y; and occupancy nj, have to be determined 
for each lattice site individually. We can diagonalize the Hamil- 
tonian (Equation 1) by the unitary transformation 


1371 


^ ^ * ^T 
Cio = F (Hino 7 OVing Y n ) (5) 


n 


where «9 are quasiparticle fermionic operators with eigenvec- 
tors Ujng and Ving. This leads to the Bogoliubov-de Gennes 
(BdG) equations 


N 

Hy Dy Sy 0 Wy 

Yint * * JT JP 
Vind Dij Ai 4 0 5j Vint 

En) yr "RS 

in} j Sij (0) Hy D; jn 

Vint N A x Vat 

0 Si Dj Ait jn 


where Dj; = ò;jU%;, and the single-particle term is given by 


E o0 
Hijo = t4, jy - (lis -078;9)8;; + S; 
with the spin-orbit coupling term 


oo’ : ^00 | ^ 
Sij =- [axs J 32a: 
l 


Here, 55^ and sge (where © is opposite to c) correspond to 
in-plane and out-of-plane spin-orbit field, respectively, and 
satisfy Sg? =(S5°)*. 


Solving numerically the BdG equations (Equation 6) we can de- 


termine the local order parameter y; and occupancy Mjo 


ti X puso e) (Ex) qn 


Nig = Ys f(&)* luus f (G5, )l (8) 


n 


where f(o) = [1 + exp(o/kgT)] !. In what follows, we shall 
inspect the spin-resolved local density of states 


Pin (0) = F| nol 8(o-6,)* vss 8(0 5, ) 


For its numerical computation we replace the Dirac delta func- 
tion with the Lorentzian function à(o) = Q[n(o? + (2)] with a 
small broadening G = 0.01 t. We have solved the BdG equations, 
considering a single magnetic impurity in a square lattice, com- 
prising Na x Np = 41 x 41 sites. We assumed U/t = —3, wt = 0, 
and determined the bound states for two representative values 


of the spin-orbit coupling X upon varying J. 
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The magnetic potential has substantial influence on the local 
order parameter yo. In particular, at some critical value J, this 
quantity discontinuously changes its magnitude and sign (see 
the upper panel in Figure 1), signaling a first-order phase transi- 
tion [28-30]. This quantum phase transition at J, is an artifact of 
the classical spin approximation. When spin fluctuations are 
allowed, a Kondo-like crossover is obtained instead of a first- 
order phase transition [31,32]. In general, the quasiparticle 
spectrum at the impurity site is characterized by two bound 


states +Eygp inside the gap A of the superconducting host 
(displayed in the bottom panel of Figure 1). These energies 


+Eysp and the related spectral weights depend on J. At J = Je 
the YSR bound states cross each other Eysr(Je) = 0 and their 
crossing signifies the ground-state parity change [33] from 
BCS-type (spinless) to the singly occupied (spinful) configura- 
tions [8,15,21,34]. Let us remark that this quantum phase transi- 
tion is also accompanied with a reversal of the YSR polariza- 
tion (see bottom panel in Figure 1). A similar behavior can be 
observed also for multiple impurities, at several critical values 
of J [35]. 


Within the BdG approach we can inspect spatial profiles of the 
YSR states by integrating the spectral weights 


x o2 
Po = f Pio (@)do 


in the interval œ € (65,05) capturing the quasiparticles at nega- 


tive/positive energies £Eysg [36]. Figure 2 illustrates the results 
obtained for à = 0 (left panel) and X. = t (right panel). We clearly 
notice a fourfold rotational symmetry (typical for the square 
lattice) and the spatial extent of YSR states reaching several 
sites away from the magnetic impurity. The non-vanishing 
difference of the spectral weight uini —|uiny? at the positive 
energy œ = +E ysp and of ving? “Win? at the negative energy 
@ = —Eysg implies the effective spin-polarization of the bound 
states (their polarization is illustrated in the bottom panel of 
Figure 1). 


For a quantitative estimation of the spatially varying magnetiza- 
tion (driven by the particle-hole asymmetry) we have com- 
puted the displaced moving average p^(r), which corresponds 
to an averaged spectral weight contained in a ring of the radius r 
and a small half-width 67. This quantity is sensitive only to the 
radial distance r from the magnetic impurity, averaging the 
angular anisotropy. Our results, presented in Figure 3, clearly 
indicate the spatial particle-hole oscillations p*(r) of the YSR 
states (compare the blue and red lines). Such particle-hole 
oscillations decay exponentially with r in agreement with 
previous studies [11,37,38]. The dominant (particle or hole) 
contributions to the YSR bound states are displayed by the 
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U(coycor) 


w/A 


por(w) — po, (o) [a-u.] 


Figure 1: The local order parameter obtained at zero temperature for weak A/t = 0.1 (red line) and strong spin—orbit coupling A/t = 1 (blue line). The 
bottom panel shows the energies and magnetic polarization po:(t») - po,(w) of YSR states obtained in the weak-coupling limit A/t = 0.1. 
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Figure 2: Spatial profiles of the YSR states y pi; obtained for |J| < Je in the absence of spin-orbit coupling (left panel) and for strong in-plane cou- 
pling À = t (right panel). The spin—orbit field is chosen along the x-axis and leads to an additional imaginary hopping term along the y-axis, which elon- 
gates the YSR states in the y-direction. The impurity spin is oriented along the (0, 0, 1) direction. 


alternating color of the background in Figure 3. We notice that 
the spin-orbit coupling seems to suppress these particle-hole 


oscillations. 


Summarizing this section, we point out that the quantum phase 
transition at J, depends on the spin—orbit coupling X and it has 
experimentally observable consequences in the magnetization 
induced near the impurity site. For weak magnetic scattering 


|J| < J, the impurity is partly screened, whereas for stronger 
couplings |J| > Je the impurity polarizes its neighborhood in the 
direction of its own magnetic moment. Similar effects have 
been previously discussed in [21], but here we additionally 
consider the role of spin—orbit coupling. First of all, such inter- 
action shifts the quantum phase transition (to larger values of J) 
and secondly it enhances the spatial extent of YSR states and 
gradually smoothes the particle-hole oscillations. 
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Figure 3: Hole-like (blue line) and electron-like (red line) displaced moving average p*(r) as a function of the radial distance r from the impurity site 
obtained for |J| < Je using dr = 0.5a. The blue and red background color indicates the dominant type (hole or particle) of the YSR states at a given dis- 
tance r. The left and right panels correspond to A = 0 and A = f, respectively. 


Magnetically polarized Majorana quasiparti- 
cles 

In this section we increase the number of impurities. Let us now 
imagine a nanoscopic chain of magnetic impurities (for instance 
Fe atoms) deposited on the surface of a conventional s-wave 
superconductor. We study the magnetically polarized bound 
states, focusing on the proximity-induced nontrivial supercon- 
ducting phase. In practice, the quasiparticle spectrum can be 
probed within STM-type setups, by attaching a conducting 
[39,40], superconducting [41], or a magnetically polarized tip 
[42]. We assume the spin-orbit interaction aligned perpendicu- 
larly to the wire and the magnetic field parallel to it, leading to 
the effective intersite pairing of identical spins and (under spe- 
cific conditions) inducing zero-energy end modes resembling 
Majorana quasiparticles. This issue has been recently studied 
very intensively but here we simply focus on the spin-polarized 
aspects of this problem. 


Due to the spin-orbit interaction, momentum and spin are no 
longer *good" quantum numbers. By solving the problem 
numerically, however, we can estimate the percentage with 
which the true quasiparticles are represented by the initial spin. 
We have recently emphasized [43], that the amplitude of inter- 
site pairing (between identical spin electrons) differs several 
times for 1 and | sectors. This leads to an obvious polarization 
of the YSR and Majorana quasiparticles (the latter appearing 
near the nanochain edges). 


Let us consider the STM-type geometry relevant to the recent 
experimental situation addressed by A. Yazdani and co-workers 
[42], which can be described by the following Hamiltonian 


^ a A prox a 
H = Htip + Hchain + Htip-chain- (9) 


We assume here that the STM tip describes a polarized fermion 
gas 


s -— 
Hn=), SkNCkoNCkoN 
k,o 


where the energy Ey = Ek — Hyo can be controlled by some 
finite detuning of the chemical potentials ux: — pw. Individual 
atoms of the nanochain are coupled with such STM tip through 


cu n Sas * ot 3 
Htip-chain = Y Vas d; kon + V; kpCkoNAi,o |- 
k,o 


For simplicity, we assume constant couplings 


Is = 2n» V, xp j &(o - &g . 
k 


The low-energy physics of such proximitized Rashba nanowire 
can be described by [44] 


A~ prox 


chain — (t = 9j Jah odio E FRashba + Hzeeman E Hipika (10) 


ijo 


where at) annihilates (creates) an electron of spin o at site i 
with energy £;, and fj; is the hopping integral. The effective 
intersite (p-wave) pairing is induced through a combined effect 
of the Rashba and the Zeeman terms 
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FRashba --—Q ` Le fio”) 


io,o 


d; s F He. (11) 


oo’ 


glpB 


Hzeeman = 


X dh(e) 4s. 


i,0,0" 


The proximity effect, which induces the on-site (trivial) pairing, 
can be modelled as [45] 


Hprox =a; (dt dt +d, 14) (13) 


with the local pairing potential A; = 'g/2. 


Figure 4 shows evolution of the spin-dependent spectrum pjg(@) 
as a function of a varying magnetic field. At a critical value 
(B = 0.2) we observe the emergence of zero-energy quasiparti- 
cles, whose spectral weights strongly depend on the spin c. 


DOS {a.u.] 


Bjt 


e 
ce 
DOS [a.u.] 


Figure 4: The effective quasiparticle spectrum pjg(w) as a function of a 
magnetic field B aligned along the nanochain obtained for o = 1 (upper 
panel) and o = | (bottom panel). The magnetic field B is expressed in 
units of t/(gup/2). 


For a better understanding of the polarized zero-energy quasi- 
particles, we present in Figure 5 the spatial profiles of the zero- 
energy (Majorana) quasiparticles. As usually such quasiparti- 
cles emerge near the edges of a nanoscopic chain, practically 
over 10 to 15 sites (see inset). Note the substantial quantitative 
difference between these zero-energy quasiparticles appearing 
in f and | spin sectors. This “intrinsic polarization” of the 
Majorana modes has been previously suggested in [46], and 
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recently we have proposed [47] their empirical detection by 
means of selective equal-spin Andreev reflection (SESAR) 
spectroscopy. 
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Figure 5: Magnetically polarized spectrum p; ;(w) - p; (w) obtained at 
w = 0 for peripheral sites of the Rashba chain. 


The main idea is to apply a bias voltage V between the STM tip 
and the superconducting substrate, inducing a charge transport 
that, in a subgap regime (|V |«A/|e]) originates from the 
Andreev (particle to hole) scattering mechanism. The polarized 
Andreev current can be expressed by the Landauer-Büttiker 


formula 


CAPS < [do T? (o)| f (o- eV) - f(@+eV)], (14) 


where transmittance is defined as 


T? (o) =r} Kid, sy +r Kid, sy 


and 
o 2 jJ 4 2 
TP (o) 2 T& dida . 


TS (9) =T} Kya - 


The anomalous Green’s functions can be computed 
numerically from the solution of the Bogoliubov-de Gennes 
equations of this model (Equation 10). The net spin current 
pm (V)= if (V) — 1, Y (V) turns out to be predominantly sensi- 
tive to the Majorana end-modes. Its differential conductance 
GP! (y) — (0/ GV I?P" (V) can thus distinguish the polarized 
Majorana quasiparticle (near V — 0) from the YSR states 
(appearing at finite voltage). 
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Bound states can leak to other side-attached nanoscopic objects. 
This proximity effect has been also predicted for the Majorana 
quasiparticles by E. Vernek et al. [48] and it has been indeed 
observed experimentally by M. T. Deng and co-authors [49]. 
Inspired by this achievement, extensive studies have been 
carried out regarding the YSR states coalescencing into the 
zero-energy Majorana states in side-coupled quantum dots 
driven by electrostatic or magnetic fields [50-52]. This issue 
would be particularly important when attempting to braid the 
Majorana end modes, e.g., in T-shape nanowires upon turning 
on and off the topological superconducting phase in its seg- 
ments. We briefly analyse here the polarized zero-energy Majo- 
rana modes leaking into the multi-site quantum dot (comprising 
ten lattice sites) side-attached to the proximitized Rashba chain 
discussed above. 


Figure 6 displays the spatial profile of the polarized spectrum 
obtained at œ = 0 as a function of the gate voltage Vg, which 
detunes the energies Vg = €; — u of the multi-site (1 < i € 10) 
quantum dot. For numerical calculations we used the model pa- 
2t, A; = 0.2t and B > Be, which guar- 
antee the Rashba chain to be in its topologically nontrivial 


rameters A = 0.155, u 


superconducting phase, hosting the zero-energy Majorana 
quasiparticles (intensive black or red regions). We clearly 
observe that for some values of Vg these Majorana modes 
spread over the entire quantum dot region. By inspecting 
Figure 6 we furthermore notice the pronounced spatial oscilla- 
tions of these zero-energy modes. In our opinion, this is a signa- 
ture of a partial delocalization of the polarized Majorana quasi- 
particles. Surprisingly, this process seems to be less efficient in 


site 
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Figure 6: Leakage of the spin-polarized Majorana quasiparticles from 
the topological superconducting phase of the Rashba chain (i 2 10) 
onto the side-attached multi-site (i € (1;10)) quantum dot. The upper 
and bottom panel show pjg(w) at w = 0 for ? and | spin, respectively. 
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the minor spin (o = |) section. This effect has to be taken into 
account, when designing nanostructures for a controllable 
spatial displacement of the Majorana modes (criticial for the re- 
alization of quantum computations with use of the Majorana- 
based qubits) either by electrostatic or magnetic means. Some 
proposals for such nanodevices have been recently discussed by 
several authors [52,53]. 


In summary of this section, we emphasize that the Majorana 
modes coalescing from the YSR states in the proximitized 
Rashba nanowire are characterized by their magnetic polariza- 
tion. Indeed, such a feature has been recently observed by STM 
spectroscopy with use of a polarized tip [42]. We have studied 
here the evolution of the polarized quasiparticle states with 
respect to the magnetic field (Figure 4) and investigated the 
spatial oscillations of the Majorana zero-energy modes near the 
chain edges (Figure 5). Finally, we analyzed leakage of the 
polarized Majorana modes on the multi-site quantum dots, 


revealing their partial delocalization (Figure 6). 


Majorana vs Kondo effect 

In previous section we have discussed the polarized Majorana 
modes leaking into side-attached objects, such as single impuri- 
ties or segments of normal nanowires. In this section we shall 
focus on the correlation effects [54-56], confronting the Majo- 
rana quasiparticle with the Kondo effect (both manifested at 
zero energy). This can be practically achieved using STM-type 
configurations sketched in Figure 7. In particular, we consider 
the subgap Kondo effect, effectively driven by the Coulomb 
repulsion U and coupling of the quantum dot (QD) with the 
normal lead Ty in presence of electron pairing (induced via I's), 
which has a significant influence on the spin-polarized bound 
states of the QD. The basic mechanism of this subgap Kondo 
effect showing up near the quantum phase transition has been 
earlier considered by us in absence of the Rashba nanowire 


[s 
Aca — An 


EN C NN 
a 


n2 


Figure 7: Schematic illustration of the quantum dot (QD) coupled be- 
tween the metallic (N) and superconducting (S) leads and hybridized 
with the Rashba nanowire, hosting the Majorana quasiparticles n4 and 
n2 at its edges. 
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[24,57]. Our considerations can be practically verified within 
STM geometry [39,40] using magnetic atoms (e.g., Fe) and 
side-coupled nonmagnetic atoms (for instance Ag or Au) 
deposited on the superconducting substrate (such as Pb or Al) 
probed with a conducting STM tip [42]. 


The topological superconducting phase, hosting the Majorana 
modes, can be driven in semiconducting wires [58,59] or in 
nanochains of magnetic atoms [39-42] through nearest-neighbor 
equal-spin pairing. The efficiency of such p-wave pairing 
differs for each spin [47], giving rise to polarization of the 
Majorana quasiparticles, with noticeable preference for the f 
sector (see Figure 4). In order to study the correlation effects we 
shall assume here a complete polarization of the Majorana 
quasiparticles. We thus focus, for simplicity, on the topological 
state originating from intersite pairing of only 1 electrons and 
consider its interplay with the correlations. Let us remark, how- 
ever, that the superconducting lead mixes both the QD spins 
with the side-attached Majorana quasiparticle [60]. In conse- 
quence we shall observe an interesting and spin-dependent rela- 
tionship between the Majorana and Kondo states that could be 
probed by the polarized Andreev (particle-to-hole conversion) 


mechanism. 


Our setup (Figure 7) can be described by the following 
Anderson-type Hamiltonian 


A= > (Hp +Hp-op )+ Hop +Hmop, (15) 


B=S,N 


where Hn corresponds to the metallic electrode, Hs refers to 
the s-wave superconducting substrate and the correlated QD is 
modeled by Hep = Y edido UR fi, where e denotes the 
energy level and U stands for the repulsive interaction between 
opposite spin electrons. The QD is coupled to both B = N,S 
reservoirs through Hg-QD = Dio pde kop +H.c.) and we 
assume a wide bandwidth limit, using the constant couplings Is. 
It can be shown [61-64] that for energies | @|«A the super-con- 
ducting electrode induces the static on-dot pairing 


Hs +Hs.op ~ Hyox = Y edid, * UR f 
o 


- (4.4, + afd). 


Taking into account the finite magnitude of superconducting 


gap [50] does not affect the main conclusions of our study. 


The effective Majorana modes of the nanowire can be modeled 
by [65] 
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ay cs AA Sum ^^t 
HMQD 7 iE mM? +A (às +4, ) 


where ù; = al are Hermitian operators and €m corresponds to an 
overlap between Majoranas. We recast these operators by 
the standard fermionic ones [66] fj =(1/ 2f fŻ) 
and ù =(-i/ 42x £4. Finally, the Hamiltonian of Equa- 


tion 15 simplifies to 


R= Fin + finon +Y edid, e Ui ie 2 (dnd, dà) 

i (16) 
^ A ^ £ 
jJ n)-m. 


es f f + tn (a, 


with the auxiliary coupling tp -A 4 A2. The subgap Kondo 
physics originates in this model from the Coulomb term Uñ, ñy 
and the effective spin-exchange interactions due to Hy.op. It 
has been shown [23,24] that under specific conditions the 
on-dot pairing can cooperate with the subgap Kondo effect. 
This particular situation occurs only near the quantum phase 


transition. 


Let us examine how the subgap Kondo effect gets along with 
the Majorana mode. Earlier studies of the correlated quantum 
dot coupled to both normal (conducting) electrodes indicated 
that the side-attached Rashba chain leads to a competition be- 
tween the Kondo and Majorana states [67-72]. For a suffi- 
ciently long wire (€m = 0) the Kondo effect persists only in the 
spin-channel |, whereas for 1 electrons there appears a dip in 
the spectral density at œ = 0. The resulting tunneling conduc- 
tance is then partly reduced (from the perfect value 2e?/h) to the 
fractional value 3e2/2h [67,68,71-73]. In contrast, for the short 
Rashba wires (with €,, # 0) the Kondo physics persists in both 


spin channels. 


In our present setup (Figure 7) the correlated quantum dot is be- 
tween the metallic and superconducting reservoirs, therefore the 
Kondo effect is additionally affected by on-dot pairing. Its in- 
fluence is mainly controlled by the ratio U/T s and partly by the 
level £, determining whether the QD ground state is in the 
spinful or spinless configuration [23,24,62,64,74]. Obviously 
the latter one cannot be screened. For instance, for the half- 
filled QD (e = —U/2) the spinful (doublet) configuration occurs 
in the regime U > I's. 


For studying the correlations we adopt perturbative treatment of 
the Coulomb potential, treating it self-consistently to the second 
order in the normal and anomalous channels [62,75]. Specific 
expressions have been provided by us in [24]. Figure 8 shows 
the spectral function po(«) for both spins obtained at zero tem- 
perature for the Coulomb potential U, covering the (spinless) 
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singlet and (spinful) doublet configurations. In the weak inter- 
action regime we observe appearance of two YSR states. For 
U zx T's these peaks merge, signaling the quantum phase transi- 
tion. The Kondo effect shows up only in the correlated limit 
(U > T's), but its spectroscopic signatures are qualitatively dif- 
ferent for each of the spins. Leakage of the Majorana quasipar- 
ticle suppresses the low-energy states of 1 electrons. We notice 
that the initial density (for tm = 0) is reduced by half, whereas 
we observe a constructive influence of the Majorana quasipar- 


ticle on opposite-spin | electrons. 


Pr 


Pi 


Figure 8: The polarized spectral function pg(w) obtained at zero tem- 
perature for the half-filled QD (€ = -U/2), Fs = 2L w, tm = 0.10 and 
several values of the Coulomb potential U (as indicated). Energies are 
expressed in units of My. 


Figure 9 shows evolution of the spectral function p;(@) for 
various couplings fm. In the weak-coupling limit we clearly 
observe a reduction (by half) of the initial density of states. 
With increasing f,, the spectrum develops the three-peak struc- 
ture that is typical for the *molecular" limit. This behavior indi- 
cates that the Majorana and Kondo states have rather a compli- 
cated relation, which is neither competitive nor cooperative. In 
fact, some novel scaling laws have been recently reported by 
several authors [69,70,76-79] also considering the correlation 
effects directly in the Rashba nanowire. 


Conclusion 
We have studied the polarized bound states of magnetic impuri- 


ties embedded in an s-wave superconducting material, taking 
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Figure 9: The spectral function p;(w) of the half-filled quantum dot 
(€ = -UI2) obtained at T = 0 for l's/TN = 2, U/TN = 4 and several values 
of tm (as indicated). 


into account the spin-orbit and/or Coulomb interactions. 
We have shown that spin-orbit coupling strongly affects the 
subgap states, both of the single impurities and their conglom- 
erate arranged into a nanoscopic chain. For the case of single 
magnetic impurity the spin-orbit interaction (i) shifts the 
quantum phase transition towards higher magnetic coupling 
Je, (ii) enhances the spatial size of the YSR states, and 
(iii) smoothes the particle-hole oscillations. For the magnetic 
chain spin—orbit coupling combined with the Zeeman term in- 
duce the topologically nontrivial superconducting state and 
indirectly give rise to substantial polarization of the Majorana 
modes (Figure 4), the oscillations of which show up near the 
chain edges (Figure 5). The polarized Majorana quasiparticles 
can also leak into other side-coupled objects, such as single or 
multiple quantum impurities (Figure 6). These polarized Majo- 
rana quasiparticles can be controlled by a magnetic field or by 
an electrostatic potential. This would be important for future 
quantum computers using qubits based on topologically pro- 
tected Majorana states. Finally, we have also confronted the 
Majorana quasiparticles with the subgap Kondo effect, 
revealing their complex relationship that can be hardly regarded 
as competitive or collaborative in some analogy to the Kondo 
effect originating from multiple degrees of freedom [80]. The 
aforementioned spin-polarized effects can be experimentally 
verified by polarized ballistic tunneling or by using STM spec- 
troscopy, relying on the selective equal-spin Andreev reflec- 
tions. 
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Abstract 


CrossMark 


We study in-gap states of a single magnetic impurity embedded in a honeycomb monolayer 
which is deposited on superconducting substrate. The intrinsic spin—orbit coupling induces the 
quantum spin Hall insulating (QSHI) phase gapped around the Fermi energy. Under such 
circumstances we consider the emergence of Shiba-like bound states driven by the 
superconducting proximity effect. We investigate their topography, spin-polarization and 
signatures of the quantum phase transition manifested by reversal of the local currents 
circulating around the magnetic impurity. These phenomena might be important for more 
exotic in-gap quasiparticles in such complex nanostructures as magnetic nanowires or islands, 
where the spin—orbit interaction along with the proximity induced electron pairing give rise to 


topological phases hosting the protected boundary modes. 


Keywords: bound states, Yu-Shiba-Rusinov states, superconductivity 
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1. Introduction 


Even a tiny content of impurities introduced to insulating 
and semiconducting materials can tremendously affect their 
charge transport, contributing particle/hole carriers from the 
donor/acceptor level to the conduction/valence band. This is in 
contrast with completely opposite (destructive) effect played 
by the magnetic impurities in superconductors where they 
break the Cooper pairs, leading to formation of the bound 
Yu-Shiba-Rusinov (YSR) or briefly Shiba states inside the 
energy gap [1]. These in-gap states can eventually activate the 
charge transport in interfaces and heterostructures, owing to 
the anomalous particle-to-hole (or hole-to-particle) Andreev 
scattering mechanism [2]. In all such cases impurities are 
intimately related with existence of the subgap states, whose 
nature differs depending on the host material. One may hence 
ask, whether there can be established any connection between 
such contrasting in-gap states of insulators and superconduc- 
tors? 

A promising platform for addressing this question could 
be a graphene sheet deposited on s-wave superconducting 
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substrate. Electrons of such carbon atoms layer reveal a num- 
ber of unique properties. Besides their Dirac-like behavior, 
stemming simply from a honeycomb geometry, the intrin- 
sic spin—orbit coupling (SOC) can induce the QSHI phase 
[3] with the spin currents circulating along its boundaries. 
Such effect has been experimentally observed in graphene 
randomly decorated with the dilute Bi? Te; nanoparticles [4] 
and in a heterostructure, consisting of a monolayer of WTe; 
placed between two layers of hexagonal boron nitride which 
has revealed topological properties up to relatively high tem- 
peratures of about 100 K [5]. Further phenomena related 
with electron pairing arise when a graphene sheet is prox- 
imitized to superconducting material [6-11]. For instance, 
graphene deposited on aluminum films acquires superconduc- 
tivity with the effective coherence length € ~ 400 nm [11], 
whereas grown on rhenium it shows high transparency of the 
interface, with the induced pairing gap A = 330 + 10 peV [8]. 
Upon introducing impurities into proximitized graphene, there 
emerge various in-gap states, manifesting either the topolog- 
ically trivial or non-trivial phases [12]. Another system for 
investigating the bound states of magnetic impurities might be 
possible in bilayer graphene, where upon twisting the carbon 
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Figure 1. Scheme of a magnetic impurity embedded in a 
honeycomb monolayer and proximitized to s-wave superconductor. 


sheets to a small ‘magic’ angle [13, 14] or tuning the interlayer 
coupling [15] the intrinsic unconventional superconducting 
phase is induced. 

Here we investigate the properties of in-gap bound states 
formed at magnetic impurity embedded into the single hon- 
eycomb two-dimensional layer and proximitized to super- 
conductor (figure 1), discussing feasible tools to unambigu- 
ously distinguish their Shiba-type character in presence of 
the QSHI phase. This problem has recently gained a great 
deal of interest, both experimentally [16—20] and theoretically 
[21-23] because similar magnetic structures, e.g. nanowires 
[24, 25] and nanostripes [26, 27] could enable realization of 
the Majorana quasiparticles. 

The spin-orbit gap in graphene is often claimed to be rather 
small, although Sichau et al [28] have estimated its magnitude 
(by means of the resistively-detected electron spin resonance) 
to be 40 eV. Under such circumstances the superconducting 
gap might be comparable to the SOC and this would be suffi- 
cient for appearance of the in-gap bound states strictly related 
with electron pairing. In what follows we perform a systematic 
study of the Shiba states, inspecting (a) their spatial extent and 
topography, (b) magnetic polarization, and (c) observable fea- 
tures of the quantum phase (0 — 7) transition manifested by 
reversal of the orbital currents circulating around the impurity 
site. 

The paper is organized as follows. In section 2 we introduce 
the microscopic model and present the method for studying 
the bound states of magnetic impurity existing in honeycomb 
sheet. Section 3 discusses influence of the insulating and super- 
conducting phases on the in-gap quasiparticles and presents 
their detailed properties. Finally, in section 4, we summarize 
the results. 


2. Model and method 


We describe the magnetic impurity embedded in a honeycomb 
sheet (figure 1) by the tight-binding Hamiltonian 


H = Íhay I Él ł Plrashba t isses (1) 


In what follows, this impurity is treated classically 
Himp =—J Cee = cle) > (2) 


where we denote the impurity site as io, and we apply the 
Kane-Mele scenario [3] for description of the itinerant elec- 


trons 
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with the nearest-neighbor hopping f, the chemical potential ji 
(which we assume to be zero unless otherwise stated), and 
the imaginary, spin-dependent, next-nearest neighbor hopping 
amplitude Aso. The latter term is responsible for inducing the 
helical edge states. The sign vi; = +1 depends on the direc- 
tion of electron hopping between the next-nearest-neighbor 
sites (+1 for clockwise and —1 for anticlockwise). The hop- 
ping terms involve the summation over (next-)nearest (((ij))) 
(ij) neighbors. Since the substrate violates the mirror inversion 
z —> —z symmetry, we also consider the Rashba spin-orbit 
interactions 


His = ie Y d Ga x di) cj». (4) 


(ij)! 


Here s^^' is the vector of the Pauli matrices, referring to spin 


i and the vector dij connects site i with its nearest neighbor 
site j. 

Finally, we assume that the honeycomb layer is proximi- 
tized to s-wave superconductor 


By = Y (Add +hc.), (5) 


1 


which induces the on-site pairing A. For computing the 
observables of interest, we perform the Bogoliubov- Valatin 
transformation 


1 
Cio = X ie In — 00. (6) 
n 


where ' denotes summation over the positive eigenvalues, and 
numerically solve the equations 


SETS = E,%, (7) 
j 


in the auxiliary (Nambu spinor) representation; = 
n n 


Qs Wis Up vr. The matrix elements read 


tit Ag 0 A^ 
. AM y ^ 0 
i= ; (8) 
0 A' hip (AI 
A 0 AD? -Gp 


where hijo = 40 ij) — (u + c Jii, )óij + giAsovjÓ(Gj)) and an 
=ide Y (sU x dy) = OE. 

Reade discussed in this paper are obtained from numerical 
diagonalization of the Hamiltonian matrix on 40 x 40 lattice 
with the periodic boundary conditions in both directions. We 
do not consider any intrinsic pairing mechanism, assuming 
that it originates solely from the proximity effect (5). Self- 
consistent treatment of electron pairing is in general important 
[29, 30], however, in the present case it would not imply any 
significant changes of the local order parameter [31]. 
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3. Subgap quasiparticles 


For a systematic analysis of the subgap quasiparticles we shall 
start by discussing the in-gap states hosted in the insulating 
(QSHI) phase and next consider their mutation caused by the 
electron pairing A. 


3.1. Impurity bound states in QSHI phase 


Let us consider the magnetic impurity in a finite-size graphene 
layer, neglecting the superconducting substrate (A = 0). 
Figure 2 shows how the intrinsic spin—orbit interaction affects 
the low-energy quasiparticles. We notice that insulating energy 
gap of the QSHI phase grows linearly upon increasing the 
Kane-Mele coupling and, around Aso = 0.21, it saturates to 
~ 1t. Concomitantly there appear two in-gap states (purple- 
dotted lines in figure 2), which are fully spin-polarized. Similar 
bound states have been previously found for a single impu- 
rity whose magnetic moment is parallel to the graphene plane 
[32]. When impurity is close enough to a perimeter of the 
sample they have been shown to hybridize with the topolog- 
ical edge states, inducing antiresonances in the transmission 
matrix. It has been also emphasized, that the bound states 
around point impurity in a two-dimensional insulator could 
distinguish between the topological and trivial phases of the 
host material [33]. 

Bottom panel in figure 2 displays the topography of the 
occupied (E < Er) bound state for two different values of Aso. 
From careful examination of the spectral weight on the lat- 
tice sites adjacent to the impurity, we can notice an oscilla- 
tory decay of the wavefunction of the bound state. Practically 
its spatial extent does not exceed 10 atomic distances, and it 
quickly vanishes for higher magnitudes of the SOC. This loss 
of spatial extent is accompanied by the simultaneous reduction 
of the spectral weight of the bound state. Closely related effects 
have been previously pointed out for the magnetic [34-38], 
non-magnetic [12, 39, 40] and both types of the scattering 
potential as well [41—43]. 


3.2. Shiba quasiparticles 


Upon coupling the honeycomb lattice to superconducting sub- 
strate, the energy gap around the Fermi energy results from a 
combined effect of the proximity induced pairing (A Æ 0) and 
the insulating phase. In general, these phenomena are known to 
be competitive as indeed manifested by suppression of the bulk 
order parameter (cj, cit) (section 3.3). From a perspective of the 
local physics (at impurity site), however, relationship between 
the QSHI and superconducting phases is much more intrigu- 
ing. By gradually increasing the pairing potential A, what can 
be achieved e.g. by reducing the external magnetic field or 
varying the temperature, we observe development of the Shiba 
quasiparticles [1] directly from in-gap quasiparticles of the 
insulating phase (figure 3). Let us remark, that direct transition 
from the insulating to superconducting phase has been theo- 
retically considered for bulk materials within the mean field 
[44] and more sophisticated many-body methods [45]. Such 
scenario could be practically realized in variety of systems, 
e.g. thin superconducting films [46], at oxide interfaces [47], 
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Figure 2. Low energy spectrum of the honeycomb lattice with the 
single magnetic impurity as a function of the Kane-Mele coupling 
Aso, assuming J = 6t and Ar = 0, A = 0. Bottom panel: 
topography of the occupied bound state for Aso = 0.1r and 
Aso = 0.2t. 
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Figure 3. Top panel: emergence of YSR states (dotted red lines) 
from in-gap quasiparticles of the QSHI phase (dotted violet lines) 
driven by the proximity induced pairing A for J = 6r, Aso = 0.11, 
H = 0. Bottom panel: same but for u = 3/3Aso. 


in organic materials [48] and possibly in the doped cuprate 
superconductors [49]. In the present context we focus on the 
subgap Shiba-like quasiparticles, which to our knowledge have 
not been considered so far. To compare our results with less 
exotic situation, we plot in the bottom panel of figure 3 the 
same situation as in the top panel, but with a value of chemi- 
cal potential which is known to close the spin—orbit gap. The 
system is then metallic and opening the superconducting gap 
results in a picture similar to traditionally understood Shiba 
states [1]. 

Let us focus in more detail on the Shiba quasiparti- 
cles. In the present case they do not obey the original for- 
mula Eysr = +A (1 — vp&(Er)J) / (1 + 7pn(Er)J) derived 
for conventional superconductors because of the vanishing 
normal density of states in graphene p,(Er) = 0 [50, 51]. 
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Figure 4. Evolution of the subgap spectrum with respect to the 
impurity potential J obtained for A = 0.2f, Ag = 0.05ż and several 
values of Ago, as indicated. The dotted red and violet lines refer to 
the YSR states and in-gap quasiparticles of insulating phase, 
respectively. 


Figure 4 displays the quasiparticle energies obtained numer- 
ically for our model as a function of the magnetic potential 
J for several values of Kane-Mele coupling Aso. The dense 
(light-blue) dots refer to a continuum, whereas the single dot- 
ted lines represent the in-gap bound states. Amongst these in- 
gap branches we can recognize the Shiba-like quasiparticles by 
their strong variation with respect to J. In particular, at some 
critical value J, they eventually cross each other, signaling a 
qualitative changeover of the ground state [52]. This quan- 
tum phase transition (QPT) manifests itself by: sign-reversal 
of the local order parameter (0 — 7t transition), abrupt onset of 
the spin polarization (section 3.3), and by qualitative changes 
(both, in magnitude and vorticity) of the locally circulating 
currents (section 3.5). 

Our analysis indicates, that Kane-Mele coupling Aso 
affects the QPT, by (i) shifting the critical coupling Jc to 
higher values (figures 4 and 6) and (ii) leading to substantial 
changes both in topography and spatial extent of the Shiba-like 
states (section 3.4). Thus the spin—orbit interaction weakens 
the efficiency of magnetic coupling J between the impurity 
and conduction electrons. Furthermore, the Shiba states no 
longer merge with a continuum even in the extremely strong 
coupling limit J — oo, in stark contrast to behavior of mag- 
netic impurities in triangular lattice of the two-dimensional 
superconductor [21] where the Kane-Mele interaction is 
absent. 


3.3. QPT 


Let us now focus on the QPT, driven by the intrinsic SOC. 
Even though variation of Aso would be rather not feasible 
experimentally, we deem that its effect can be instructive for 
understanding mutual relationship between the on-site pair- 
ing and the spin-—orbit interaction. Bottom panel of figure 5 
presents the eigenenergies, corresponding to the same set of 
model parameters as in figure 2 but in presence of finite A 
and Ag. We observe that energy gap of superconducting states 
(~ 0.20) gradually evolves into the gap of QSHI which sat- 
urates around Aso œ 0.2t. We have selected strong enough 
magnetic coupling J = 6t which allows for the QPT driven by 
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Figure 5. The polarization P and order parameter at the impurity 
site and for the bulk as functions of Aso (top panel). QPT driven at 
Aso œ 0.05t (bottom panel). Results obtained for the model 
parameters A = 0.2f, Ax = 0.05r, and J = 6t. 
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Figure 6. Variation of the QPT point (corresponding to crossing of 
the Shiba states) versus the Kane—Mele coupling Aso and impurity 
potential J. 


Aso. The upper panel of figure 5 displays the bulk polarization, 
defined as 


P= 59> (m) - (ru), (9) 


where nig = >>, [lu f(En. T) + |, f(—E,.T)] is the 
average number of electrons with spin c at site i, the order 
parameter at the impurity site, and its value averaged over the 
entire sample. At Aso ~ 0.051 the order parameter at impurity 
site reverses its sign and its absolute value gradually increases 
upon increasing the Kane- Mele coupling. Simultaneously the 
bulk magnetization is abruptly quenched as the system shifts 
to the unpolarized ground state. These characteristic features 
of the QPT [1] in the present case originate from the intrin- 
sic SOC. On the other hand, the bulk order parameter does 
not undergo any dramatic changes (it slowly diminishes upon 
increasing Aso). Such conditions should be taken into account 
when considering formation of the Majorana bound states at 
the ends of magnetic chains deposited on the proximitized 
honeycomb sheet [25]. 

The shift of Je with increasing Aso is displayed as a phase 
diagram in figure 6. The black continuous line denotes the 
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Figure 7. Spatial distribution of the occupied (negative value) Shiba 
quasiparticle obtained for A = 0.2t, Ar = 0.05t, J = 4t, using 

Aso = 0 (left panel) and Aso = O0.1* (right panel). The density of 
states p;(Esniba) is normalized with respect to the largest value in the 
bottom panel. 


critical coupling Je at different values of Aso. Initially the 
shift of QPT is not meaningful, but starting from Aso — 
0.03t we observe onset of a linear variation. This increase 
also points out, that the spin—orbit interaction suppresses the 
effective coupling of the impurity spin with the conduction 
electrons [36]. 


3.4. Topography of Shiba quasiparticles 


Let us now inspect the real-space shape (topography) of 
the Shiba states. Figure 7 presents spatial maps of the 
density of states at the energy of electron-like (occupied) 
bound state, both in absence and in presence of the intrin- 
sic spin-orbit interaction. One can see that without the 
Kane—Mele interaction, the topography of Shiba state has its 
usual character with exponential variation of the wavefunction 
~ exp —r/8. 

We remark, that spectral weight is differently distributed 
in each sublattice. Close to the impurity site ry = (0,0) most 
of the spectral weight of the Shiba quasiparticles appears in 
the B-sublattice sites, whereas further away the A-sublattice 
(in which the impurity resides) gains more and more spec- 
tral weight. Also the rotational symmetry of the topographic 
shape reveals a bipartite character. Close to the impurity site 
the shape has a C3 rotational symmetry, reflecting the fact 
that each site has three nearest-neighbors (cf bottom panels 
in figure 2), whereas at larger distances, the spectral weight 
distributed in the A sublattice resembles a hexagon with Cg 
rotational symmetry. Precise evaluation of the bound states 
wavefuntions in this case would be a challenging task for 
future experimental measurements. Topography of the bound 
states changes dramatically, when the intrinsic SOC is taken 
into account. Bottom panel in figure 7 illustrates a strong ten- 
dency towards localization of the Shiba states in vicinity of 
the magnetic impurity. Their spectral weight is spread over a 
few adjacent sites and we no longer observe any preference 
for dominance of only one sublattice. These properties of the 
Shiba states resemble the features typical for in-gap quasiparti- 
cles of magnetic impurity embedded in a non-superconducting 
QSHI. Such reduction of the spatial extent could be impor- 
tant for engineering the topologically non-trivial phases, as e.g. 
chain of magnetic impurities can host the Majorana quasiparti- 
cles only when the bound states of dilute impurities hybridize 
to form a Shiba band capable of undergoing the topological 
phase transition. 


0.4 1.0 
10 

0.3 0.8 

0.6 
Qy 02 = 
E 0.4 a ) 

-10 94 — 10 0.2 =1, 

0.0 0.0 = 


S Glodzik and T Domański 


3.5. Local currents 


Another signature of the QPT in our system can be seen by cur- 
rents induced around the magnetic impurity [53]. We compute 
" 


the local charge flow, using the Heisenberg equation i 55 = 


(Ini, H)). Setting the convention h = 1, and ignoring the ae 
which merely induce on-site fluctuations of charge, we obtain 


_ ii 
= E («d Ch Cje) — (cow) 


To > (vis? 7 (Cha'Cio) — vas?" (c},sCio)) 
oa'((j)) (10) 


BARI |(se" x dj). ch cjo) 
oo" (j) 
— [uu x di). ( cas) . 


Applying the Bogoliubov—Valatin transformation (6), and 
T: 
making use of the fact that if ®; = (ui Wi Ups Vj) is the 
eigenvector of matrix (8) with an eigenvalue E,, then ®; = 
T : : 
(— UE UD.—u a is also an eigenvector of the same 
matrix, but with an eigenvalue —E,, we get 


= Bs uu Uj, frp(En) = c.c. ) 
jon 


T Y Ox" Uig i, fro(En) + c-c.) an 
(j)oo'n 
+ Aso D ise” ugu ujo fro(En) + c. ehi 
((j))oo'n 


where A^ defined in section 2. 

Figure 8 shows the real-space maps of microscopic currents 
and the maximum magnitude of bond current in the system 
with respect to the impurity coupling strength J. We empha- 
size, that reversal of these currents vorticity (compare the top 
panels of figure 8) occurs exclusively when the system is in 
the non-trivial QSHI phase. Explanation of this behavior could 
be the following. It has been shown in reference [36] that the 
intrinsic SOC supresses the effective coupling J of impurity 
with the conduction electrons. We have observed that with 
Aso = 0 the sites belonging to the same sublattice as the impu- 
rity site polarize easily in the direction of the magnetic moment 
of the impurity. This effect is pronounced only for J > Jc, as 
more sites around the impurity align their magnetic moments. 
The situation changes with increasing SOC which weakens 
the effective impurity coupling. For small J, the magnetic 
momentis hardly screened by the closest neighboring sites and 
becomes more efficient when the coupling exceeds the critical 
value Jc, forcing the neighboring sites to align their magneti- 
zation with the impurity. This in turn reverses the direction of 
the current. The Shiba states (labeled YSR in red vector plots 
in figure 8) are the ones that cross the Fermi energy during the 
QPT, hence only their contribution to the current shows this 
change of direction, in contrast to the bound states (BS in blue 
vector plots in figure 8) discussed in section 3.1. Those states 
hardly change their energy with increasing J, and their contri- 
bution to the current does not change during the QPT. Bottom 


J. Phys.: Condens. Matter 32 (2020) 235501 


[Jt-25] Bs Jft — 85 
-y —-N 
AEN » LN * 
i a Si d ar bi 
EN T 
CUESS] vss LESE: 
` i os eges 3 
/ \ "2 \ H 
y "uit y ise e 
ES CR ET .* 


mm 


Imaz 


, | 
-10. — Je -5 0 5 
| Jjt 


Figure 8. Vector maps of the currents around the magnetic impurity 
obtained for J = 5t (left) and J = 8.5t (right) presenting 
contributions of QSHI bound states (BS) and Shiba states (YSR). 
Bottom panel shows the maximum current versus the coupling J 
when taking into account the whole spectrum (black dashed line), 
only QSHI bound states and Shiba states (blue and red lines 
respectively). Other parameters: A = 0.2t, Aso = 0.11, 

Ar = 0.050. 


panel in figure 8 presents the maximum value of the current in 
the system. When summing over every state n in equation (11), 
the current drops discontinuously at QPT. This is because as 
can be observed from the contribution of only the Shiba states 
(red) and the QSHI bound states (blue), after the reversal of 
current direction of the Shiba states, both contributions com- 
pete, and the effective maximum current is greatly reduced. 
Detection of such orbital patterns might be performed using 
an integrated quantum imaging platform where graphene sheet 
is connected to an array of the atomic-sized magnetic sensors 
[54, 55] or local conductivity atomic force microscopy suit- 
able for probing electronic current paths with a diameter in the 
nanometer range [56]. 


4. Conclusions 


We have theoretically investigated the energetic, magnetic and 
topographic features of in-gap quasiparticles of a single mag- 
netic impurity embedded in the graphene sheet and prox- 
imitized to the superconducting substrate. We have shown 
that subtle interplay between the intrinsic spin—orbit inter- 
action (responsible for the energy gap of the QSHI phase) 
and the proximity-induced electron pairing enables emergence 
of the Shiba-type quasiparticles directly from in-gap states 
of the insulating phase. We have discussed in detail this 
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intriguing behavior and proposed several methods for its 
empirical verification. 

Furthermore we have found, that upon varying either the 
magnetic coupling J (feasible in STM experiments [18]) or 
strength of the spin—orbit coupling Aso a pair of the Shiba 
bound states could cross at the Fermi energy, causing quantum 
phase transition of the ground state. This usually leads to sign- 
change of the local order parameter [1], but in the present situ- 
ation it would be also uniquely manifested by a reversal of the 
vorticity and abrupt change of the absolute magnitude of the 
local currents circulating around the impurity site. Our numer- 
ical calculations have additionally revealed, that the spin—orbit 
coupling pushes such QPT crossing towards much higher val- 
ues of J and substantially reduces the extent of Shiba states 
(similar to the in-gap states of insulating phase), partly affect- 
ing their topographic patterns. We have carefully inspected 
their spatial profiles in each sublattice of the graphene sheet. 

We hope that phenomena discussed here for the single-site 
magnetic defects [57, 58] might stimulate further consider- 
ations of the topological insulating and/or superconducting 
phases in more complex magnetic structures, like nanowires 
[24, 25], nanoscopic islands [19] or stripes [26, 27], where the 
Majorana-type quasiparticles can be realized. It would be also 
worth to extend our study on the quantum impurities, address- 
ing the subgap Kondo physics of the conventional [59, 60] and 
topological [61] superconductors. 
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Abstract 

Recently it was discovered that superconductivity in transition metal dichalcogenides (TMDs) is 
strongly affected by an out-of-plane spin—orbit coupling. In addition, new techniques of fabricating 
2d ferromagnets on van der Waals materials are rapidly emerging. Combining these breakthroughs, 
we propose a realization of nodal topological superconductivity in TMDs by fabricating nanostruc- 
tured ferromagnets with an in-plane magnetization on the top surface. The proposed procedure does 
not require application of external magnetic fields and applies to monolayer and multilayer (bulk) 
systems. The signatures of the topological phase include Majorana flat bandsthat can be directly 
observed by scanning tunneling microscopy techniques. We concentrate on NbSe; and argue that the 
proposed structures demonstrating the nodal topological phase can be realized within existing 
technology. 


1. Introduction 


According to the modern approach to condensed matter physics, novel states of matter can be realized in 
designer systems by combining simpler building blocks. This view implies that our access to new phases of 
matter and emergent quantum particles is ultimately limited only by our imagination and ability to manipulate 
matter. The designer approach to 1d topological superconductivity [1, 2] has stirred enormous activity, aiming 
to fabricate Majorana quasiparticles [3, 4] and harness them in applications. While most of the previous work 
has targeted gapped phases, here we propose fabrication ofa 2d nodal phase with flat Majorana edge bands. 

Our proposal is based on two recent breakthroughs, the observation of Ising superconductivity in transition 
metal dichalcogenides (TMDs) [5-12] and the discovery of novel techniques to fabricate 2d magnetic structures on 
van der Waals materials [13]. Due to the inversion-breaking structure of monolayers, the quasiparticles experience 
strong valley-dependent out-of-plane spin—orbit coupling (SOC). As a result, superconductivity exhibits 
remarkable robustness in the presence oflarge magnetic fields inducing a Zeeman splitting far exceeding the Pauli 
limit. Theimportance ofthe SOC in bulk 2H layer structures has been long overlooked probably because the bulk 
has inversion symmetry. While the stacked monolayers that make up the bulk exhibit staggered SOC that restores 
inversion symmetry as a whole, thelayers are weakly coupled and quasiparticles in individual layers are subject to 
strong Ising SOC [14]. This is particularly interesting since the pioneering work by Kane and Mele [15] identified 
an Ising type SOC in graphene asa crucial ingredient ofthe quantum spin-Hall phase. 

The rich electronic properties of TMDs have stimulated several proposals for topological and 
unconventional superconductivity [16—26]. In [27, 28] it was proposed that an Ising spin—orbit coupled 
monolayer TMD (e.g. NbSe;) in the presence of an in-plane magnetic field could realize a topological 
superconducting state characterized by a nodal bulk and flat Majorana edge bands. In our work we consider how 
the state can be observed in magnetic structures fabricated on top of an Ising superconductor, treating NbSe, asa 
particular realization. We show how magnetic nanostructures on 2H-NbSe, give rise to the nodal topological 
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phase signaled by the appearance of flat Majorana edge bands. In the light ofthe ground-breaking success in 
fabricating 2d magnetic nanostructures on van der Waals materials, our design takes advantage of the latest 
advances and has several crucial advantages over the magnetic field induced phase [29-32]. Our setup does not 
only remove the necessity of external fields but, importantly, relaxes the requirement of manufacturing 
monolayer or few layer systems. The proposed topological state engineering works equally well for bulk systems 
since disruptive orbital effects arising from the in-plane magnetic fields are completely absent in our design. In 
multilayer systems the topological state is formed in the surface layer in the area in contact with the magnetic 
structure. The further advantages ofthe proposed setup include the possibility of fabricating well-defined 
nanostructures of topological elements. This comes with the benefit that the Majorana flat bands can be directly 
observed by scanning tunneling microscopy (STM) by studying the local density of states (LDOS) on magnetic 
islands. The STM measurements could be employed in resolving the spatial structure ofthe flat bands, thus 
providing a smoking-gun signature ofthe nodal topological phase. Considering that the process of growing 
magnetic islands on top of NbSe; systems is already experimentally on the way, we expect that our predictions 
will find experimental confirmation in the near future. 


2. Model 


Asa particular model ofa TMD we consider a tight-binding representation of 2H-NbSe». This layered structure 
consists of stacked units of Nb atoms on a triangular lattice sandwiched by Se layers. To study the appearance and the 
properties ofthe nodal topological phase, it is convenient to device a minimum phenomenological model that 
faithfully displays the essential features ofa real material. Since the relevant bands near the Fermi energy mostly 
derive from Nb d orbitals, the band structure can be qualitatively reproduced by a model with one orbital per site on 
atriangular lattice. This approximation treats the system effectively as quasi-2d structure, neglecting the less 
important Se-derived 3d bands near the T point. While the tight-binding parameters and the band structure vary 
with the number of layers, the quasi-2d bands remain largely unchanged due to the weak interlayer coupling. Also, 
experimental observations of magnetic impurities in bulk 2H-NbSe; highlight the 2d nature of the magnetic subgap 
spectrum and the qualitative validity of treating the layers as effectively decoupled [33]. Thus, while our model is 
expected to be most faithful in monolayer systems, the experimental evidence suggests it can reasonably capture the 
qualitative behavior ofthe top layer in contact with the magnet in multilayer systems. In a recent experiment 
examining Shiba states of Fe impurities in NbSe2 ([34]), it was observed that the subgap energy of single magnetic 
impurity states may get modulated by the presence ofthe charge density wave (CDW). In the magnetic island 
structure we study (see figure 1 for a schematic view), similar effect could perhaps give rise to modulation of 
magnetically induced subgap bands with CDW periodicity. However, this is not expected to be crucial for the flat 
edge bands. The Hamiltonian on a triangular lattice with the zigzag edge parallel to the x direction can be written as 
Ê = Ain + Hsoc + Hy + Asc, where 


i m t t t 
Hy, — —t y» CiCjg m y n ty CiCjo» 
(jho ic (ijh 
Asoc = -i D eio cid, 
v 
Ay =— > M, (io? Gus 
Lag 
Ago = Y Alici + H.c.). (1) 
i 


Symbols t and t; correspond to the nearest and next-nearest neighbor hopping and p is the chemical 
potential. Additionally, \ parametrizes the out-of-plane Ising-type SOC, the magnetic material gives rise to 
magnetic splitting M, and A describes the superconducting pairing. Pauli matrices c act on the spin degrees of 
freedom. The symbols (i, j) and ( (i, j)) denote the summation over nearest and next-nearest neighbors 
respectively, 6j; is the angle the vector connecting i and j sites makes with the positive x axis (so that e? = +1). 

The system is most conveniently analyzed by passing to momentum space and working in the Nambu basis 
W = (e cp Gr cl p where the Bogoliubov-de Gennes Hamiltonian becomes 


H = Eo(k)7; + Eso(k)o; + Mo, + Atyoy. (2) 


The additional set of Pauli matrices 7 act in the particle-hole space. The normal and spin-orbit hopping terms 
are given by 


E(k) = -2| cos (k,a) + 2cos (= ) cos Ga ) an] cos (ky 43a) + 2cos (s ) cos Es ) 
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Figure 1. Schematic view ofthe studied system. A magnetic island with in-plane polarization (blue circle) on a top surface of NbSe; 
structure. Triangular lattice sites indicate the positions of Nb atoms. The in-plane TMD structure determines the positions ofthe flat 
bands. This is independent on the magnetization direction. 


M 


Figure 2. Momentum space characteristics ofa monolayer TMD. (a) Fermi surface and the Brillouin zone with high symmetry 
directions as dashed purple lines. Parameters used: t; = t, jj = —0.8t, A = 0.1tand b = 43a, where ais the lattice constant. (b) The 
evolution of point nodes along a high symmetry I'-M line in the superconducting state with the values of Zeeman field in inset boxes. 
Same parameters as (a) but with A = 0.1t. The inset in the last panel shows a close-up of the purple dashed box. 


and 
Eso(k) = 20 sin (kya) sin(5* St) sin (5 za (3) 
Diagonalization of the Hamiltonian reveals four energy bands: 
E? = Ed + Eg, + M? + A +2 EED + M2) + Mê. (4) 


In the next section we show how the essential normal state features of NbSe, as well as the nodal topological 
superconducting state emerge from this model. 


3. Nodal topological superconductivity and flat edge bands 


3.1. Normal state properties 

To model the system in the normal state without magnetization, we first set M = A = Oand plot the Fermi 
surface in figure 2, for a set of parameters that qualitatively reproduces the numerically calculated and 
experimentally resolved band-structure of NbSe; [14, 35]. We observe the spin-orbit split [ pocket and more 
visibly split pockets around K and K’ points. Within the second-neighbor hopping model that we use, we are 
also able to resolve a slight trigonal warping of the pockets. The SOC polarizes the electron spins at Kand K’ 
points in the opposite out-of-plane directions. The splitting vanishes along l'-M that define mirror symmetric 
lines. Note that in a multilayer system the spin—orbit coupling is staggered and A has opposite sign in adjacent 
layers. Nevertheless, the top layer relevant for the topological nodal state in multilayer systems will experience a 
strong Ising SOC just as a monolayer system. 


3.2. Nodal topological state 

In the presence of superconductivity and magnetization the spectrum is given by equation (4). The spectrum is 
symmetric with respect to zero energy and branches with the minus sign in front of the square root will 
determine the properties near Fermi energy. For A > M the system is fully gapped while for M > A the system 
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Figure 3. Results in the infinite strip geometry. Top panel: E(k,) spectrum of the nanoribbon. The number of atoms parallel to the 
zigzag direction N = 4001. The red dashed line marks the range of k, for which we expect a flat band to form. Bottom panel: winding 
number equation (5). The inset shows a close-up of the red dashed box. Parameters used: t; = t, p 0.85, A = 0.14, A = 0.14 
M, = 2A. 

y 


becomes gapless with isolated nodal points ko satisfying E(ko) = 0. From the dispersion (4) we can calculate that 
at M = Aa pair of nodes (with opposite winding) nucleate on the crossing point of the Fermi surface and the 
I-M lines. This will give rise to six pairs of nodes in the Brillouin zone. By increasing M the nodes of opposite 
winding move away from the Fermi surface along the T'—M lines. This evolution ofthe nodes is presented in 
figure 2(b). The nodes cannot be gapped out by small lattice-symmetric perturbations, hinting to a topological 
nature ofthe phase. The topological character will be made more explicit below when we discuss the appearance 
ofthe flat edge bands. 

In analogy to the Fermi arc surface states that connect the surface projections of bulk band-touching points 
in topological semimetals, the 2d nodal phase supports edge states that connect the edge projections ofthe 
nodes. Due to superconductivity, the edge states in the studied system have a Majorana character. The dispersion 
ofthe edge bands is flat and the bands appear on edge terminations for which the projections ofthe nodes with 
opposite winding do not cancel. Atomic positions on a single layer of NbSe; form a hexagonal lattice with 
triangular Nb and Se sublattices. The flat band is the most prominent for an armchair edge termination while it 
shrinks to a point for a zigzag edge. In order to probe the flat bands, we envision the studied system in an infinite 
ribbon geometry. We assume periodic boundary conditions (PBCs) in y direction, yielding a flat band with 
maximal extension. So obtained strip spectrum E(k,) is presented in the top panel of figure 3. The large 
mismatch ofthe relevant energy scales A ~ 0.01t would require considering system sizes that are impractical for 
numerical calculations. Therefore we employ exaggerated parameters A ^ 0.1tfor real-space calculations. 
However, below we argue that the flat bands remain well-defined even for realistic parameters (see figure 6). 
Close to the nanoribbon BZ edge we see a perfectly flat band at zero energy which connects the nodes. Closer to 
kyb/v ~ +0.5, there exists another flat band connecting opposite nodes. For the employed parameters, the gap 
between the bulk states and the edge band is much smaller and the states are less localized to the edge. As a result, 
the flat band exhibits oscillating departures from the gap center. As shown below, the direct evaluation ofthe 
topological invariant protecting the flat band confirms that the inner edge bands exhibit perfect flatness in the 
large system limit. 

To elucidate the topological nature of the nodal phase, we calculate the topological invariant protecting the 
flat bands. This will also provide a definite connection between the edge bands and the edge projection ofthe 
nodes. The Hamiltonian (2) anticommutes with matrix C = T, dx, hence it belongs to the class BDI, and can be 
characterized by the winding number in odd spatial dimensions. In a strip geometry we can Fourier transform 
the 2d tight-binding Hamiltonian in y direction. Thus, the resulting partially transformed Hamiltonian 
describes hopping in a 1d chain perpendicular to the strip while k, is regarded as a parameter. For some intervals 
of k, the 1d Hamiltonian is topologically non-trivial and the 1d chain hosts end states. The flat band can be 
regarded as the union ofthe end states ofthe Id Hamiltonian. The winding number as a function of k, can be 
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Figure 4. Flat bands are formed between edge projections of nodes with opposite chirality. The edge parallel to y direction supports 
four edge bands (two shown) with maximal extension. The band that connects the projections oftwo pairs of nodes is doubly 
degenerate. On the edge parallel to x axis the opposite nodes project on top of each other and the flat band vanishes. 
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Figure 5. The local density of states w.r. to the energy E and the site index in the nanoribbon geometry. Parameters same as in figure 3, 
expect for A = 0.3t. 


obtained by evaluating the invariant [36] 
1 27 
v(k) = — f dk, Tr[, a H- 10, H], (5) 
271 J—20 


and is illustrated for the studied model in the bottom panel of figure 3. We can see that close to the ribbon BZ 
edge, where there is a visible flat band, the value of the invariant is (kj) = 1. There is also an interval of non- 
trivial momenta around k,b/m = 0.5, with higher value of the topological invariant 1(k,) = —2. The different 
signs and values of the winding number come from the addition of nodes when projecting them onto the k, 
direction as depicted in figure 4. Physically this means that the flat band is doubly degenerate when |v (k,)| = 2. 
Close to the BZ edge, one pair of nodes adds up and the invariant is equal to unity. The mid-gap bands close to 
the middle of the BZ that suffer from the proximity of the bulk bands, coincide with invariant value 1(k,) = —2. 
Thus they are unambiguously identified as flat bands in the large system limit. 

To see that the flat bands are indeed localized at the edges of the nanoribbon, we calculated the LDOS as a 
function of energy and the site index in the direction perpendicular to the PBC 


A, E) = 5 [ting (ky) 76 Œ = En (ky) F Wino (ky) P6CE sp E,(k,)], (6) 


kyn,o 


where we sum over all values of momentum, every state n and both spin directions c = },|.The Dirac functions 
are approximated by Lorentzians with broadening 0.0011. The result presented in figure 5 shows a strong 
enhancement ofthe zero-energy LDOS at the edges ofthe strip. 


3.3. Nodal phase on magnetic islands 

The most direct link between our theory and experiments is probably provided by TMD systems with magnetic 
islands grown on top. These systems, expected to display rich interplay of competing orders, are currently 
becoming accessible in experiments. The study of coexistence of superconductivity and magnetization has along 
history. A single magnetic impurity induces a pair of bound states in the superconducting energy gap. Their 
shape in real space reflects the symmetry ofthe underlying lattice, and can extend far away from the impurity 
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Figure 6. Evolution ofthe point nodes (crossing points) and low-lying spectrum for different values of A/t. (a) ‘Inner’ nodes 
(b) ‘outer’, close to the Brillouin zone edge nodes. Dotted lines in (b) denote the position of flat bands in a finite system. While the 
length ofthe flat band is diminished, the gap to the bulk states remains significant. 


Figure 7. Local density of states integrated over the two lowest-lying states, expressed as the size ofthe yellow dots. Dashed lines 
outline the shapes of magnetic islands. Parameters used: t; = t, ji 0.8, à = 0.15, A = 0.35 M, = 2A. 


depending on the dimensionality of the substrate and other factors like e.g. spin—orbit interaction [33, 37]. In the 
case of a collection of impurities, a Shiba band forms which can undergo a topological phase transition in the 
deep-dilute limit. Such two-dimensional islands of magnetic impurities are recently receiving much attention in 
extensive theoretical and experimental studies [38—41], with proposals of engineering hybrids with different 
dimensionality [42]. In this context, the present work turns a new leaf on topological state engineering through 
magnetization, generalizing the previous efforts to obtain gapless topological phases. 

For real space calculations we have constructed a triangular lattice consisting of 85 x 86 atomic sites and 
applied PBCs in both directions. Then we have applied Zeeman field pointing in the y direction on a collection of 
sites that comprise the magnetic island. Our approach offully diagonalizing the lattice model is laden with high 
computational cost, hence the limited finite size ofthe lattice. Due to the finite size limitations, the flat band 
physics can be properly resolved only for the two lowest-lying (closest to zero energy) states, and for the toy 
parameters, in which the superconducting energy gap is strongly exaggerated (A ~ t). However, this is sufficient 
to demonstrate qualitatively the effects of the nodal phase. To demonstrate that we show in figure 6 how the low- 
lying spectrum between the nodal points evolve as values of A/tare decreasing. One can see that the distance 
from the gap center to the bulk bands is approximately the same for different values of ^, and for the outer nodes 
its value remains a significant fraction of the superconducting gap. This means that the flat bands are always 
clearly resolved from the bulk states and that the calculations with exaggerated A/t values do not qualitatively 
change the results. The summed LDOS of the low-energy states is shown in figure 7. In panel (a) we assumed a 
circular shape ofthe magnetic island with radius r = 28a from the center of the lattice, whereas in (b) the island 
has a rectangular shape of the size 60a x 40a, where ais the lattice constant, assumed to be equal to unity. The 
number of atomic sites comprising the island region is approximately equal in both geometries. The underlying 
lattice is periodic in every direction, hence the edge of the island is the only edge in our system. The nodal 
topological phase emerges only when there is an in-plane magnetization, so the edge of the island is the 
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boundary between a gapped trivial phase and the nodal topological phase. Because ofthe structure of point 
nodes (as explained in section 3.2) and their mutual annulment in the k, projection, the edge states will never 
appear in parallel to x direction. We observe that indeed, regardless of the shape of the magnetic island, localized 
Majorana modes appear only on the edges parallel to y direction. As stressed above, a quantitative study of the 
edge bands in real space comes with substantial computational cost and is beyond the scope of the present work. 


3.4. Note on the Rashba effect 

In general, monolayer systems on a substrate and surfaces of a bulk system break the mirror symmetry 
perpendicular to layers. Therefore a Rashba-type spin—orbit coupling may be non-zero. The Rashba effect, 
giving rise to an in-plane spin-orbit field, has adversarial effect on the nodal state as it generally gaps out the 
nodes. The Rashba contribution to equation (2) can be implemented through the nearest-neighbor hopping 
term 


Hg = En, (k) ox + En, (K)720y, ds 
with 
Ep, (k) = anv sin (55) cos (t). 


Ep, (k) = -orfino + sin (= ) cos (553)). 


where ag is the Rashba constant. As can be easily verified, the chiral symmetry protecting the edge band is 
broken anda strict topological protection is not realized. Unlike time-reversal symmetry (which is always 
present in the absence of magnetism) and particle-hole symmetry (which is exact for superconducting systems), 
chiral symmetry typically requires some degree of tuning of external conditions. A typical example is the 
direction of magnetic field. For some field orientations the system exhibits chiral symmetry while for other this 
does not hold. Since some degree of disorder is almost unavoidable in solid state systems, chiral symmetry is 
typically only approximate. However, the magnitude of breaking ofthis symmetry is important in assessing how 
detrimental it is for the observability ofthe edge modes. A small symmetry-breaking perturbation pushes the 
edge modes away from the zero energy and they acquire a weak dispersion, but they are still present. In contrast 
to the Ising SOC which is determined by the lattice structure of TMDs, the Rashba SOC is case specific. 
Therefore it can also be very weak, especially in multilayer systems where ripples do not play a role. 


4. Discussion and outlook 


The physical realization ofthe proposed system could be a multilayer or monolayer TMD in contact with 
magnetic insulating material grown on top. The material should support in-plane magnetization and not 
perturb the system significantly. A 2d magnetic insulator with a high-quality contact to TMD would be ideal for 
this purpose. In fact, the requirement to bean insulator is inconsequential since proximity effect can make a thin 
magnet superconducting. The recent breakthroughs in fabricating 2d magnets down to monolayer thickness on 
van der Waals systems provide a promising avenue for our proposal [13, 29—32]. An interesting candidate for the 
ferromagnet is VSe; which is a TMD itself and can be epitaxially grown on NbSe; and other systems of interest. 
Furthermore, it has been observed that structures based on VSe; layers exhibit an in-plane magnetization on 
different substrates [31]. Considering the emerging nature and rapid development of the field of 2d magnets, 
increasing number of material candidates are likely to emerge soon. In practice the magnetic material also 
induces a non-magnetic potential which could shift the chemical potential of substrate and thus change the 
Fermi surface. However, the precise filling is not crucial for realization of the nodal phase. It seems plausible that 
in the recent experiment [43] the magnetism in the system was not sufficiently homogeneous or that it is 
significant only near the edges ofthe VSe islands. While the first experimental result does not seem to exhibit 
the signatures ofthe flat band, it is quite possible that further development in the sample fabrication and 
experimental methods could yield the desired result. 

The nodal phase can be experimentally identified by observing the Majorana flat bands. As discussed above, 
the flat bands give rise to enhanced zero-energy LDOS on certain edge terminations on the island. An ideal probe 
toaccess this information is STM. In principle, the surface LDOS can be directly measured as a function of 
energy. This would resolve the flat bands in space and energy. As long as the in-plane magnetization is 
sufficiently strong to drive islands to the nodal phase, the flat band is most pronounced and suppressed in the 
same spatial directions for all islands irrespective of the direction of their magnetization. Observation of an 
enhanced zero energy LDOS on edges with common tangent for multiple islands would constitute a smoking 
evidence on the nodal phase. 
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In the present work we employed NbSe; as a candidate material for the topological state engineering. 
However, for the existence of the nodal phase the crucial features are the (quasi-)2d nature ofthe system, the 
Ising SOC within a layer, superconductivity and in-plane magnetization. Thus we expect that other TMDs 
would also provide promising candidates for the proposed system. 

In summary, we proposed transition metal dichalcogenides with ferromagnetic structure on top as 
promising candidates to realize nodal topological superconductivity and flat Majorana edge bands. The systems 
could be fabricated and analyzed within existing experimental techniques. 


Acknowledgments 


This work was supported by National Science Centre (NCN, Poland) grant 2017/27/N/ST3/01762 (SG). We 
would like to thank Peter Liljeroth for initial stimulus of the project and for the numerous discussions regarding 
experimental realizations. We also acknowledge fruitful discussions with Tadeusz Domański. 


ORCID iDs 


Szczepan Glodzik © https: //orcid.org/0000-0001-7928-9959 


References 


[1] Oreg Y, Refael Gand von Oppen F 2010 Helical liquids and Majorana bound states in quantum wires Phys. Rev. Lett. 105 177002 

[2] Lutchyn R M, Sau J D and Sarma S D 2010 Majorana fermions and a topological phase transition in semiconductor-superconductor 
heterostructures Phys. Rev. Lett. 105 077001 

[3] Mourik V, Zuo K, Frolov S M, Plissard S R, Bakkers E P A M and Kouwenhoven L P 2012 Signatures of Majorana fermions in hybrid 
superconductor-semiconductor nanowire devices Science 336 1003—7 

[4] Das A, Ronen Y, Most Y, Oreg Y, Heiblum M and Shtrikman H 2012 Zero-bias peaks and splitting in an Al-InAs nanowire topological 
superconductor as a signature of Majorana fermions Nat. Phys. 8 887 

[5] Saito Y etal 2015 Superconductivity protected by spin-valley locking in ion-gated MoS, Nat. Phys. 12 144 

[6] XiX, Wang Z, Zhao W, Park J-H, Law KT, Berger H, Forró L, Shan J and Mak K F 2015 Ising pairing in superconducting NbSe; atomic 
layers Nat. Phys. 12 139 

[7] XiX, Berger H, Forró L, Shan J and Mak K F 2016 Gate tuning of electronic phase transitions in two-dimensional NbSe, atomic layers 
Phys. Rev. Lett. 117 106801 

[8] LuJ M, Zheliuk O, Leermakers I, Yuan N F Q, Zeitler U, Law K T and Ye J T 2015 Evidence for two-dimensional Ising 
superconductivity in gated MoS, Science 350 1353-7 

[9] Navarro-Moratalla E et al 2016 Enhanced superconductivity in atomically thin TaS; Nat. Commun. 7 11043 

[10] Huang C etal 2018 Inducing strong superconductivity in WTe; by a proximity effect ACS Nano 12 7185-96 

[11] dela Barrera S C et al 2018 Tuning Ising superconductivity with layer and spin—orbit coupling in two-dimensional transition-metal 

dichalcogenides Nat. Commun. 9 1427 

[12] Sohn E etal 2018 An unusual continuous paramagnetic-limited superconducting phase transition in 2d NbSe2 Nat. Mater. 17 504—8 

[13] Shabbir B, Nadeem M, Dai Z, Fuhrer M S, Xue Q-K, Wang X and Bao Q 2018 Long range intrinsic ferromagnetism in two dimensional 

materials and dissipationless future technologies Appl. Phys. Rev. 5 041105 

[14] Bawden L et al 2016 Spin-valley locking in the normal state ofa transition-metal dichalcogenide superconductor Nat. Commun. 7 

HZ 

[15] Kane C Land Mele EJ 2005 Quantum spin Hall effect in graphene Phys. Rev. Lett. 95 226801 

[16] Zhou B T, Yuan N FQ, Jiang H-L and Law K T 2016 Ising superconductivity and Majorana fermions in transition-metal 

dichalcogenides Phys. Rev. B 93 180501 

[17] Shaffer D, Kang J, Burnell F J and Fernandes R M 2019 Crystalline nodal topological superconductivity in monolayer NbSe; 

arXiv:1904.05470 

[18] Yuan N FQ, MakK F and Law K T 2014 Possible topological superconducting phases of MoS, Phys. Rev. Lett. 113 097001 

[19] Hsu Y-T, Vaezi A, Fischer MH and Kim E-A 2017 Topological superconductivity in monolayer transition metal dichalcogenides Nat. 

Commun. 8 14985 

[20] Fischer MH, Sigrist M and Agterberg D F 2018 Superconductivity without inversion and time-reversal symmetries Phys. Rev. Lett. 121 

157003 

[21] Taniguchi K, Matsumoto A, Shimotani H and Takagi H 2012 Electric-field-induced superconductivity at 9.4 K in a layered transition 

metal disulphide MoS, Appl. Phys. Lett. 101 042603 

[22] Oiwa R, Yanagi Y and Kusunose H 2018 Theory of superconductivity in hole-doped monolayer MoS, Phys. Rev. B 98 064509 

[23] Móckli D and Khodas M 2018 Robust parity-mixed superconductivity in disordered monolayer transition metal dichalcogenides Phys. 

Rev. B98 144518 

[24] Sosenko E, Zhang J and Aji V 2017 Unconventional superconductivity and anomalous response in hole-doped transition metal 

dichalcogenides Phys. Rev. B 95 144508 

[25] Liu C-X 2017 Unconventional superconductivity in bilayer transition metal dichalcogenides Phys. Rev. Lett. 118 087001 

[26] Chen W and An J 2019 Phys. Rev. B 100 054503 

[27] He W-Y, Zhou B T, HeJ J, Yuan N FQ, Zhang T and Law K T 2018 Magnetic field driven nodal topological superconductivity in 

monolayer transition metal dichalcogenides Commun. Phys. 1 40 

[28] WangL, Rosdahl T O and Sticlet D 2018 Platform for nodal topological superconductors in monolayer molybdenum dichalcogenides 

Phys. Rev. B98 205411 


IOP Publishing 


New J. Phys. 22 (2020) 013022 S Glodzik and T Ojanen 


[29] O'Hara D J etal 2018 Room temperature intrinsic ferromagnetism in epitaxial manganese selenide films in the monolayer limit Nano 
Lett. 18 3125-31 

[30] Gong C etal 2017 Discovery of intrinsic ferromagnetism in two-dimensional van der Waals crystals Nature 546 265 

[31] Bonilla M, Kolekar S, Ma Y, Diaz H C, Kalappattil V, Das R, Eggers T, Gutierrez H R, Phan M-H and Batzill M 2018 Strong room- 
temperature ferromagnetism in VSe; monolayers on van der Waals substrates Nat. Nanotechnol. 13 289-93 

[32] Huang B et al 2017 Layer-dependent ferromagnetism in a van der Waals crystal down to the monolayer limit Nature 546 270 

[33] Ménard G C et al 2015 Coherent long-range magnetic bound states in a superconductor Nat. Phys. 11 1013 

[34] Liebhaber E, Gonzalez S A, Baba R, Reecht G, Heinrich B W, Rohlf S, Rossnagel K, von Oppen F and Franke K J 2019 Yu-Shiba- 
Rusinov states in the charge-density modulated superconductor NbSe; arXiv:1903.09663 [condmat.mes-hall] 

[35] Wijayaratne K, Zhao J, Malliakas C, Chung D Y, Kanatzidis bc M G and Chatterjee U 2017 Spectroscopic signature of moment- 
dependent electron-phonon coupling in 2H-TaS, J. Mater. Chem. C 5 11310 

[36] Volovik G E 2009 The Universe in a Helium Droplet (International Series of Monographs on Physics) (Oxford: Oxford University Press) 
[37] Ptok A, Glodzik S and Domański T 2017 Yu-Shiba-Rusinov states of impurities in a triangular lattice of NbSe; with spin—orbit 
coupling Phys. Rev. B 96 184425 

[38] Róntynen] and Ojanen T 2015 Topological superconductivity and high Chern numbers in 2d ferromagnetic Shiba lattices Phys. Rev. 
Lett. 114 236803 

[39] Róntynen J and Ojanen T 2016 Chern mosaic: topology of chiral superconductivity on ferromagnetic adatom lattices Phys. Rev. B 93 
094521 

[40] Palacio-Morales A, Mascot E, Cocklin S, Kim H, Rachel S, Morr D K and Wiesendanger R 2019 Atomic-scale interface engineering of 
Majorana edge modes in a 2d magnet-superconductor hybrid system Sci. Adv. 5 eaav6600 

[41] Ménard G C, Guissart S, Brun C, Leriche R T, Trif M, Debontridder F, Demaille D, Roditchev D, Simon P and Cren T 2017 Two- 
dimensional topological superconductivity in Pb/Co/Si(111) Nat. Commun. 8 2040 

[42] Kobialka A, Domanski T and Ptok A 2019 Delocalisation of Majorana quasiparticles in plaquette-nanowire hybrid system Sci. Rep. 9 
12933 

[43] Kezilebieke S et al 2019 Electronic and magnetic characterization of epitaxial VSe; monolayers on superconducting NbSe; 
arXiv:1909.10208 [cond-mat.mtrlsci] 


arX1v:2002.02141v2 [cond-mat.mes-hall] 8 Feb 2020 


Topological superconductivity in a designer 


ferromagnet-superconductor van der Waals heterostructure 


Shawulienu Kezilebieke, * Md Nurul Huda,! Viliam Vaiio,! Markus 
Aapro,! Somesh C. Ganguli,! Orlando J. Silveira,! Szczepan Glodzik,? 


Adam S. Foster,’ Teemu Ojanen,*^?* and Peter Liljeroth! * 


! Department of Applied Physics, Aalto University School of Science, 
PO Box 15100, 00076 Aalto, Finland 
? Institute of Physics, M. Curie-Sktodowska University, 20-031 Lublin, Poland 
34 WPI Nano Life Science Institute (WPI-NanoL SI), 
Kanazawa University, Kakuma-machi, Kanazawa 920-1192, Japan 
4Computational Physics Laboratory, Physics Unit, 
Faculty of Engineering and Natural, Sciences, 
Tampere University, PO Box 692, FI-33014 Tampere, Finland 
? Helsinki Institute of Physics PO Box 64, FI-00014, Finland 
(Dated: February 11, 2020) 


Abstract 

'The designer approach has become a new paradigm in accessing novel quantum phases of mat- 
ter. Moreover, the realization of exotic states such as topological insulators, superconductors 
and quantum spin liquids often poses challenging or even contradictory demands for any single 
material. For example, it is presently unclear if topological superconductivity, which has been 
suggested as a key ingredient for topological quantum computing, exists at all in any naturally oc- 
curring material. T'his problem can be circumvented by using designer heterostructures combining 
different materials, where the desired physics emerges from the engineered interactions between 
the different components. Here, we employ the designer approach to demonstrate two major 
breakthroughs — the fabrication of van der Waals (vdW) heterostructures combining 2D ferromag- 
netism with superconductivity and the observation of 2D topological superconductivity. We use 
molecular-beam epitaxy (MBE) to grow two-dimensional islands of ferromagnetic chromium tri- 
bromide (CrBr3) on superconducting niobium diselenide (NbSe2) and demonstrate the existence of 
the one-dimensional Majorana edge modes using low-temperature scanning tunneling microscopy 
(STM) and spectroscopy (STS). The fabricated two-dimensional vdW heterostructure provides a 


high-quality controllable platform for electronic devices harnessing topological superconductivity. 


There has been a surge of interest in designer materials that would realize electronic 


L7. Topological superconductors are 


responses not found in naturally occurring materials 
one of the main targets of these efforts and they are currently attracting intense attention 
due to their potential as building blocks for Majorana-based qubits for topological quantum 
computation’ ?. Majorana zero-energy modes (MZM) have been reported in several different 
experimental platforms, with the most prominent examples being semiconductor nanowires 
with strong spin-orbit coupling and ferromagnetic atomic chains proximitized with an s-wave 


7,9-14 


superconductor . It is also possible to realize MZMs in vortex cores on a proximitized 


15-17 or on FeTeo.55Seo.45 superconductor surface!*:?, In these 


topological insulator surface 
cases the MZM were spectroscopically identified as zero energy conductance signals that 
are localized at the ends of the one dimensional (1D) chain or in the vortex core. In two- 
dimensional systems, 1D dispersive chiral Majorana fermions are expected to localize near 
the edge of the system (Fig. 1A). For example, it was proposed that the dispersing Majorana 
states can be created at the edges of an island of magnetic adatoms on the surface of an s- 
wave superconductor”? ??, Experimentally, promising signatures of such 1D chiral Majorana 
modes have recently been reported around nanoscale magnetic islands either buried below 
a single atomic layer of Pb??, or adsorbed on a Re substrate?4, and in domain walls in 
FeTeg559eo9457. However, these types of systems can be sensitive to disorder and may 
require interface engineering through, e.g., the use of an atomically thin separation layer. In 
addition, it is difficult to incorporate these materials into device structures. These problems 
can be circumvented in van der Waals (vdW) heterostructures, where the different layers 
interact only through vdW forces!. VdW heterostructures naturally allow for very high 
quality interfaces and a multitude of practical devices have been demonstrated. While vdW 
materials with a wide range of properties have been discovered, ferromagnetism has been 
notably absent until recent discoveries of atomically thin CraGesTeg?9, CrI3?7 and CrBrs?5?9, 
The first reports relied on mechanical exfoliation for the sample preparation, but CrBr3°° 
and Fe3GeTes?! have also been grown using molecular-beam epitaxy (MBE) in ultra-high 


vacuum (UHV). This is essential for realizing clean edges and interfaces. 


Among the various known vdW materials, the recently discovered monolayer ferromagnet 
transition metal trihalides combined with transition metal dichalcogenide (TMD) supercon- 
ductors form an ideal platform for realizing 2D topological superconductivity (Fig. 1A). 


Here, we use MBE to grow high-quality monolayer ferromagnet CrBr3 on a NbSes super- 
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FIG. 1. Realization of topological superconductivity in CrBra3-NbSe» heterostructures. 
(A) Schematic of the experimental setup. (B,C) Schematic of the bandstructure engineering 
to realize topological superconductivity. Effect of adding spin-orbit interactions and weaker and 
stronger Zeeman-type magnetization on the low-energy band structure in the normal (B) and 
superconducting states (C). (D) STM image of a monolayer thick CrBrs island grown on NbSe» 
using MBE (STM feedback parameters: Vpias = +1 V, J = 10 pA, scale bar: 10 nm). (E) Atomically 
resolved image on the CrBr3 layer (STM feedback parameters: Vpiag = +1.7 V, J = 0.5 nA, image 
size: 19 x 19 nm?). (F) Calculated structure and the induced spin-polarization from density- 
functional theory calculations. (G) Experimental d//dV spectroscopy on the NbSe2 substrate 
(blue), the middle of the CrBrg island (red) and at the edge of the CrBrs island (green) measured 


at T = 350 mK. 


conducting substrate. The mirror symmetry is broken at the interface between the different 
materials and this lifts the spin degeneracy due to the Rashba effect. Therefore, we have 
all the necessary ingredients — magnetism, superconductivity and Rashba spin-orbit cou- 
pling — required to realize a designer topological superconductor????. We demonstrate the 
existence of the one-dimensional Majorana edge modes using low-temperature scanning tun- 
neling microscopy (STM) and spectroscopy (STS). Realizing topological superconductivity 
in a van der Waals heterostructure has significant advantages compared to the other possi- 
ble platforms: vdW heterostructures can potentially be manufactured by simple mechanical 
exfoliation, the interfaces are naturally very uniform and of high quality, and the structures 
can be straightforwardly integrated in device structures. Finally, layered heterostructures 
can be readily accessed by a large variety of external stimuli making external control of 2D 
topological superconductivity potentially possible by electrical?*^, mechanical?, chemical®®, 
and optical approaches?’. 

Pioneering theoretical works???? demonstrated that topological superconductivity may 
arise from a combination of out-of-plane ferromagnetism, superconductivity and Rashba- 
type spin-orbit coupling, as illustrated in Fig. 1B,C. In this scheme, the Rashba coupling 
lifts the spin-degeneracy of the conduction band while Zeeman splitting due to proximity 
magnetization lifts the remaining Kramers degeneracy. Adding superconductivity creates a 
particle-hole symmetric band structure and the superconducting pairing opens gaps at the 
Fermi energy. In our theoretical model for magnetically covered NbSeg, a similar picture 
arises for the real band structure around any of the high symmetry points of the hexagonal 
Brillouin zone (T, K, or M) where Rashba coupling vanishes. Depending on the magnitude 
of the magnetization induced gap M and the position of the Fermi energy u, the system 
enters a topological phase when |e(ko) — u| X M, where e(ko) is the energy of the band 
crossing at the high symmetry point in the absence of magnetization. This is due to the 
created effective p-wave pairing symmetry. The power of the designer approach with vdW 
heterostructures comes from the fact that the different components retain their intrinsic 
properties allowing for rational design of emergent quantum phases of matter and realizing 
schemes such as shown in Fig. 1B,C. 

In Fig. 1D, we show a constant-current STM image of the CrBrs island grown on a freshly 
cleaved bulk NbSes substrate by MBE (details of sample growth and STM experiments are 
given in the Supplementary Material (SM)). The CrBrs islands show a well-ordered moiré 
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superstructure with 6.3 nm periodicity arising from the lattice mismatch between the CrBrs 
and the NbSes layers. Fig. 1E shows an atomically resolved STM image of the CrBr3 mono- 
layer, revealing periodically spaced triangular protrusions. These features are formed by 
the three neighboring Br atoms as highlighted in the Fig. 1F (red triangle) showing the 
fully relaxed geometry of CrBr3/NbSes heterostructure obtained through density functional 
theory (DFT) calculations (see SM for details). The measured in-plane lattice constant is 
6.5 A, consistent with the recent experimental value (6.3 A) of monolayer CrBr3 grown on 
graphite? and our DFT calculations. As expected for a weakly-interacting vdW heterostruc- 
ture, DFT calculations further confirm that the CrBr3 monolayer retains its ferromagnetic 
ordering with a magnetocrystalline anisotropy favouring an out-of-plane spin orientation as 
shown in Fig. 1F. The magnetization density (Fig. 1F) shows that the magnetism arises 
from the partially filled d orbitals of the Cr+ ion. While the largest magnetization den- 
sity is found close to the Cr atoms, as expected, there is also significant proximity induced 
magnetization on the Nb atoms in the underlying NbSe; layer. 

This establishes that our system our system has the required key ingredients for topolog- 
ical superconductivity and now we probe the resulting emergent quantum matter with STS 
measurements at a temperature of T = 350 mK. Fig. 1G shows experimental d/ /dV spectra 
(raw data) taken at different locations indicated in Fig. 1D (marked by filled circles). The 
dI/dV spectrum of bare NbSes has a hard gap with an extended region of zero differential 
conductance around zero bias, which can be fitted by a double gap s-wave BCS-type spec- 
trum (see SM for details). In contrast, the spectra taken in the middle of the CrBrs island 
have small but distinctly non-zero differential conductance inside the gap of the NbSes sub- 
strate. Moreover, we observe pairs of conductance onsets at 40.3 mV around zero bias (red 
arrows). The magnetization causes the formation of energy bands (dubbed Shiba bands) 


that exist inside the superconducting gap of the substrate???, 


By introducing spin-orbit 
interactions (as discussed above), the system can be driven into a topological phase with 
associated closing and reopening of the gap between the Shiba bands. 

We observe edge modes consistent with the expected Majorana modes along the edge of 
the magnetic island that are the hallmark of 2D topological superconductivity????4, More- 
over, the spectroscopic feature of the Majorana edge mode appears inside the gap defined 


by the Shiba bands (the topological gap) and is centred around the Fermi level (Ep). A 
typical spectrum taken at the edge of the CrBr3 island is shown in Fig. 1G, where a peak 
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FIG. 2. Electronic structure of CrBr3-NbSe2 heterostructures. (A) The band structure 
of the spin-split Nb d-band used in the effective model for topological superconductivity with 
magnetization M — 100 meV. The inset shows the 1** Brillouin zone, where the six M-points 
and the Rashba texture around them has been highlighted. (B) Calculated phase diagram of the 
magnetized NbSes based on the effective low-energy model. The color scale indicates the energy 
gap Egap (in the units of A). (C) The calculated topological gap A; as a function of the Rashba and 
magnetization energies (in units of the superconducting gap A). (D) Calculated band structure of 


the topological phase based on a phenomenological tight-binding model (see SM for details). 


localized at Ep is clearly seen. Furthermore, there are two side peaks located at +0.41 mV 
(see SM for a more detailed analysis) that are very close in energy to the Shiba bands on 
the CrBrs layer. 

To account for the experimental observations and to corroborate the topological nature 
of the edge modes, we model the system through DFT calculations and develop an effective 
low-energy model for the system (see SM for details). The band structure of the Nb d-states 
derived band used in the effective model is shown in Fig. 2A (direct comparison with DFT 
is shown in the SM). Topological superconductivity can be generated when magnetization is 


sufficiently strong to push one of the spin-degenerate bands at a high-symmetry point above 
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the Fermi energy. We identify the observed topological phase as a state arising from the 
gap-closing transition at M point with a Chern number C — 3. The calculated topological 
phase diagram in Fig. 2B indicates that for a reasonable magnetization of M < 100 meV, 
the C = 3 state lies approximately 100-200 meV below the Fermi energy of pristine NbSes. 
Thus, a small renormalization of the chemical potential resulting from the contact with the 
magnetic material will drive the system to the topological phase. The two other nontrivial 
phases that originate from gap closings at the I point and K points give rise to topological 
phases with C — —1 and C — —2. Realization of either of these phases would require notably 
larger shifts in chemical potential (~ 0.6 eV), making them improbable for the experimental 
observations. The absolute values of the nontrivial Chern numbers can be understood by a 
three-fold rotational symmetry (see SM). 

The key quantity characterizing robustness of the nontrivial phase is the topological 
energy gap A;. This scale should be much larger than temperature for the state to be 
observable in experiment. In the simple parabolic band model this quantity can be esti- 
mated by A; = akp/|(akr)? + M?]'?, where a is the Rashba coupling and kp the Fermi 
wavelength?". The calculated gap based on our more realistic tight-binding (TB) model 
is shown in Fig. 2C. Based on the experimental results shown in Fig. 1G, the topological 
gap is A, £ 0.3A. The calculated band structure in a strip geometry corresponding to the 
experimental gap is shown in Fig. 2D, where we see the Majorana edge modes crossing the 
topological gap. The edge modes are seen to coexist with the bulk states in a finite subgap 
energy window in agreement with experimental observations. 

To further elucidate the properties of the Majorana edge modes, we have carried out 
spatially resolved d//dV spectroscopy over the edge of the CrBrs island (Fig. 3A). It can 
be seen that the edge mode is more localized at the edge at zero bias (deeper within the 
topological gap), but becomes more delocalized at energies closer to the topological gap edge. 
Moreover, the energy dependence of the main feature of the edge mode LDOS is such that it 
splits off from the top edge of the topological gap inside the CrBrs island, smoothly crosses 
the topological gap and merges with its lower edge outside the CrBr; island. In order to 
visualize the evolution of the spatial extension of the Majorana edge modes, we have recorded 
grid d//dV spectroscopy maps (Fig. 3B-G). At Ep, the Majorana edge modes are confined 
within ~ 2.4 nm of the edge of the island (see SM). In addition to the edge mode signature 


close to the Fermi level, there is also enhanced LDOS at the energies similar to size of the 
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FIG. 3. Spatially resolved spectroscopy of the Majorana zero modes. (A) dI/dV spec- 
troscopy over the edge of the CrBr3 island (STM topography shown on the top). (B-G) STM to- 
pography and spatially resolved LDOS maps extracted from grid spectroscopy experiments. STM 
feedback parameters: (A) Vbias = +1 V, I = 10 pA; (B) Vias = +0.8 V, I = 10 pA. Scale bars: 
(A) 4 nm; (B) 12 nm. (H-N) Corresponding calculated linespectra (H) and LDOS maps (J-N) 


with the island shape shown in (I) (see SM for details). 


topological gap (Fig. 3E,F) where we also see significant excitations inside the magnetic 
island. The theoretically computed LDOS (Fig. 3H-N, see SM for details) reproduces the 
essential features of the experimental results and shows exponential localization of Majorana 


edge modes at the edge of the island at subgap energies. In particular, we confirm the 
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experimental finding that the distribution of the spectral weight of the edge mode along 
the edge (excluding the regular moiré pattern) should be non-uniform. This stems from 
the geometric irregularities of the island boundary with characteristic length scale that 
is comparable to the edge mode penetration depth. However, this does not imply the 
edge modes are discontinuous along the edge. It simply means that the interference effects 
near edge irregularities suppress the visibility of the edge mode due to finite experimental 
resolution. The simulations also reproduce the slightly dispersing peaked LDOS across the 
island edge, and the enhanced LDOS at the islands edges at and above the topological gap 
edge. 

It is conceivable that the experimentally observed edge modes could possess a topologi- 
cally trivial origin, unrelated to the existence of a topological superconducting state. How- 
ever, in addition to the near quantitative match with the theoretical results incorporating 
the main ingredients of the experimental system, the edge mode signature is experimentally 
very robust. We consistently observe it in our hybrid vdW heterostructures on all CrBrs 
islands, irrespective of their specific size and shape (see SM for more examples). To prove 
the observed edge modes of the hybrid heterostructures are strongly linked to the supercon- 
ductivity of the NbSes substrate, we have carried out experiments in magnetic fields up to 
4 T, suppressing superconductivity in the NbSes substrate. All features associated with the 
gap at the center of the island and the edge modes disappear in the absence of superconduc- 
tivity in NbSes (see SM for details). This rules out trivial edge modes as the cause of the 
observed results. Another non-topological reason for resonances close to the Fermi energy, 
the Kondo effect, should also be present in the normal state and can hence be ruled out as 


well. 


In conclusion, our work constitutes two breakthroughs in designer quantum materials. 
By fabricating vdW heterostructures with 2D ferromagnet epitaxially coupled to super- 
conducting NbSe2, we obtained a near ideal designer structure exhibiting two competing 
electronic orders. The induced magnetization and spin-orbit coupling renders the supercon- 
ductor topologically nontrivial, supporting Majorana edge channels which we characterized 
by STM and STS measurements. The demonstrated heterostructure provides a high-quality 
platform for electrical devices employing topological superconductivity. This is an essen- 
tial step towards practical devices employing topological superconductivity and the demon- 


strated system would in principle allow electrical control of the topological phase through 


10 


electrostatic tuning of the chemical potential. 
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We analyze the topological superconductivity, and investigate the spectroscopic properties mani- 
fested by zero-energy modes, induced in a metallic strip embedded into a Josephson-type junction. 
Focusing on the Majorana polarization of such quasiparticles we propose feasible means for its em- 
pirical detection, using the spin selective Andreev reflection method. Our study reveals a gradual 
development of a transverse gradient of the Majorana polarization across the metallic strip upon 
increasing its width. We also inspect the spatial profile and polarization of the Majorana quasipar- 
ticles in the presence of a strong electrostatic defect. We show that, depending on its position, such 
a defect can lead to a substantial localization of the Majorana mode. 


I. INTRODUCTION 


'Topological materials, including those which are either 
insulators or superconductors, differ qualitatively from 
their ordinary counterparts due to the emergence of in- 
gap modes. Such quasiparticles develop at boundaries 
or internal defects and are topologically protected (thus 
being good candidates for stable qubits), and obey frac- 
tional statistics (which is appealing for quantum compu- 
tations). Experimental efforts for the realization of these 
exotic quasiparticles have so far largely focused on one- 
dimensional structures, e.g. semiconducting nanowires 
proximitized to superconductors,'? nanochains of mag- 
netic atoms deposited on superconducting substrates,° |! 
and lithographically fabricated nanostructures.!? An- 
other direction in pursuit of topological superconductiv- 
ity relies on two-dimensional systems, where the in-gap 
quasiparticles are chiral modes.!? !* Such Majorana edge 
modes have indeed been observed in STM measurements, 
using nanoscopic islands of magnetic atoms deposited 
on superconducting substrates.!? ?? Further interesting 
perspectives are related with mixed-dimensionality sys- 
tems, where the localized and delocalized Majorana 
quasiparticles coexist with one another.?!?? [n partic- 
ular, nanowires attached to larger structures?’ could 
enable a controllable transfer of the Majorana modes 
between these constituents,?* indirectly probing their 
Chern numbers.?? 


Yet another promising platform for the realization 
of topological superconductivity hosting localized Majo- 
rana modes has been suggested in Refs.?9?" using metal- 
lic strips with strong spin-orbit coupling embedded be- 
tween two superconducting leads with differing phases 
(see Fig. 1). Signatures of zero-energy modes have been 
observed in such heterostructures, consisting both of alu- 
minium on indium arsenide?? and an HgTe quantum 
well coupled to thin-film aluminium.?? The major virtue 
of a Josephson-type geometry is its tunability to the 
topologically non-trivial regime, easily controlled exper- 
imentally by varying the phase difference. Another way 
for a controllable transition to the topological phase is 
possible by using two gate-tunable Josephson junctions 


FIG. 1. A schematic view of a metallic strip (dark purple), 
embedded between superconducting regions (yellow) which 
differ in phase by $, probed by a polarized STM tip (light 
gray). A magnetic field B= Boĉ is applied to the whole 
structure. 


(i.e. a SQUID geometry), as recently reported for epitax- 
ial Al/InAs structures.?? 

Experiments on these Josephson junction heterostruc- 
tures??? have triggered further intensive studies.?! ?* 
'The proximitized metallic strips are hoped, for instance, 
to enable a current-controlled braiding of the Majorana 
modes.?? It has also been suggested"? that weak disorder 
promotes localization of the Majorana quasiparticles. 

Intrigued by this prediction, we study here the spatial 
profiles and polarizations of the Majorana modes. We 
also consider a single strong point-like electrostatic scat- 
tering potential placed in various regions of the proximi- 
tized metallic strip. For this purpose we perform numer- 
ical calculations within the Bogoliubov-de Gennes treat- 
ment. Our study reveals that, when this local defect is 
placed in an interior of the metallic stripe its influence on 
the Majorana modes is naturally practically negligible, 
but when placed near a region of the existing Majorana 
quasiparticle we observe a tendency towards reducing the 
spatial extent of the zero-energy modes, in some analogy 
to what has been predicted by Haim and Stern.?9 These 
phenomena are in stark contrast to the properties of Ma- 
jorana quasiparticles of strictly 1-dimensional systems, 
where sufficiently strong local defects can produce addi- 
tional pairs of the Majorana modes. 

We also inspect the Majorana polarization introduced 


in Refs. 21, 37-39 and considered by one of us.4° We show 
that this quantity could be particularly useful in charac- 
terizing the zero-energy quasiparticles of metallic strips 
whose experimentally reported length-to-width ratio was 
about 20 (Ref. 28) or 100 (Ref. 29). It is known that quasi 
2-dimensionality induces transverse gradients in the Ma- 
jorana polarization.?.?9:4 Similarly quasi 1-dimensional 
systems, such as those considered here, can have long lo- 
calization length scales for the Majorana bound states 
which can also induce transverse gradients in the Majo- 
rana polarization. Strictly the gradient is induced in the 
phase of the Majorana polarization, and in the magnitude 
of the Majorana polarization relative to its density. We 
prove that the magnitude (absolute value) can be probed 
by selective equal spin Andreev reflection spectroscopy.*? 
A similar method has previously been applied for inspect- 
ing the spin polarization of the Majorana quasiparticles 
of Fe atom chains using a magnetic STM tip.!? The ad- 
vantage of the method proposed here is that it would be 
a direct probe of the Majorana nature of the quasiparti- 
cles. Transverse gradients of the Majorana polarization, 
relative to the density, are an indication of a delocaliza- 
tion process that ultimately could be detrimental to the 
zero-energy modes. 

This paper is organized as follows. In Sec. II we present 
the microscopic model and briefly outline methodological 
details. In Sec. III we investigate the spatial profiles and 
polarization of the Majorana modes, focusing on their 
evolution with respect to varying the width of the homo- 
geneous metallic strip. Sec. V discusses the Majorana lo- 
calization driven by a point-like electrostatic defect, and 
Sec. VI summarizes the main results. 


II. MICROSCOPIC MODEL 


For a description of the planar Josephson heterostruc- 
ture, see Fig. 1, we employ the microscopic scenario dis- 
cussed in Refs. 26, 27, and 43. The model Hamiltonian, 


H=Hot+tHz+Hs, (1) 
consists first of the free term 


Ho = X [iN (dij x 855). — too] d dj: — UY di, dig 

(3) ic 

0,0 

(2) 

describing itinerant electrons hopping all over the sample. 
t is the hopping integral between the nearest neighbor 
sites on a square lattice, À is the strength of the Rashba 
spin-orbit coupling, d;; is the vector connecting nearest 
neighbors, and o stands for the vector of the Pauli ma- 
trices. The second (Zeeman) term 


Hz — Bod > di ,o% odio! (3) 


accounts for the influence of an external magnetic field Bo 
which is parallel to the interface between the metallic and 


superconducting regions, as reported experimentally.?5:2? 


The last part appearing in the model Hamiltonian (1) 
describes the on-site pairing in the left (Sr) and right 
(Sg) superconducting regions, 


Hs - (dl, dl, + Hc.) , (4) 
where 
Ae7i#/2 for i € Sr, 


Ae'?/2 fori € Sp, and (5) 
0 foie N. 


A = 


Here the metallic strip region is denoted by N. The phase 
difference between the superconducting layers Sp and Sr, 
is ¢ and A is real. 

We studied the finite-size version of this model, con- 
sisting of Nz sites along x-direction and N, sites along 
y-axis. For specific computations we assumed N, = 91 
and N, = 30, unless stated otherwise. The eigenstates 
and eigenenergies of the heterostructure were determined 
numerically, solving the Bogoliubov de-Gennes equations 
with the canonical transformation 


G)-EDURSIO) @ 


n 


where XP stand for the Bogoliubov quasiparticles which 


diagonalize the Hamiltonian: H = >, Enyi Yn + const. 
In particular we have calculated the local density of 
states 


pi(w) = 5 5 (luk, ilw — En) + t ?6(w + En)) (7) 


n,o 


As we consider a finite size system we have broadened 
the delta function peaks to Lorentzian functions of width 
0.024. 

Another quantity of interest is the Majorana polariza- 
tion??? for an eigenstate |Yn), 


Pin = (bs |Cf; vn) = 5 Tog 2W; Vio , (8) 


where f; is projection onto site z and C is the particle-hole 
operator. This quantity allows one to probe directly the 
Majorana nature of the eigenstates, and its experimental 
measurement is discussed in Sec. III B. For convenience 
we introduce 


Pi = Pap Pi (9) 
where 
Pio = 2u; U2; (10) 


for the zero-energy (n = no) quasiparticles. More gen- 
erally one may wish to consider the particle-hole overlap 
Ugo UJo,: In particular the equal-spin pairing (o1 = c3) 
induced between the neighboring sites ? and j for the 
zero-energy quasiparticles (Ej, = 0 = Em) is of interest. 


More details on this issue are discussed in Sec. III B. 


III. TOPOGRAPHY OF MAJORANA MODES 


Upon substituting the metallic strip between the su- 
perconducting reservoirs, their Cooper pairs leak into the 
normal region, inducing on-site electron pairing. This 
proximity effect is efficient nearby the bulk superconduc- 
tors, up to distances smaller than the coherence length £. 
Here we consider metallic samples comprising a few Nw 
atomic rows, whose spatial width Nua < £, where a is the 
inter-atomic distance. Under such a condition the prox- 
imity effect induces superconductivity across the entire 
metallic region. The appearance of the topological su- 
perconducting phase, however, can be realized only with 
triplet pairing which can be achieved by combining con- 
ventional superconductivity with the spin-orbit Rashba 
interaction and Zeeman splitting.” 


It has been demonstrated?9?7" that a transition from 
the topologically trivial to the nontrivial superconduct- 
ing state is sensitive to the Josephson phase $. Charac- 
teristic features of the emerging Majorana quasiparticles 
can, however, additionally depend on the width N,,a of 
the metallic strip. In what follows we analyze such qual- 
itative changes and propose a method for their empirical 
detection (Sec. IV). 
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FIG. 2. The spatial profiles, p;(0), of the Majorana quasi- 
particles appearing in a metallic strip consisting of 1, 2 and 
5 rows of atomic chains, as indicated. We have used the 
model parameters A = 0.25t, à = m, A = 0.5t, Bo = 0.1t, 
u = —3.75t. The Majorana polarization (10) for the area in 
the green square is shown in Fig. 3. 


A. Zero-energy spectral function 


Focusing on the optimal condition for the topological 
superconducting state ¢ = m, we have checked that the 
on-site pairing (di, d;+) spreads nearly uniformly onto the 
metallic strip, both along and across it. Furthermore, we 
also noticed some feedback of the metallic sector onto 
superconducting regions manifested by partial reduction 
of the local pairing (sometimes referred to as the inverse 
superconducting proximity effect). Next, when combined 
with the spin-orbit interaction and the Zeeman field, the 
proximitized metallic stripe develops inter-site pairing of 
identical spin electrons, i.e. triplet pairing. For suffi- 
ciently strong magnetic fields Bo, such a triplet supercon- 
ducting phase becomes topologically nontrivial, leading 
to the emergence of the zero-energy quasiparticles.?9:27 


Fig. 2 displays the spatial profiles of the local density 
of states at zero energy p;(0), obtained for very narrow 
metallic strips. We note that the Majorana quasiparti- 
cles of the narrow metallic strip are well localized at its 
ends. Their overall topography is practically identical 
with all features of one-dimensional systems, including 
the characteristic oscillations along the metallic strip.*4 
It comes as perhaps some surprise that this narrow width 
of metallic region is neither essential for the development 
of the topological superconducting phase, nor important 
for the spatial profile of the Majorana modes. Even in 
the extreme case Nw = 0, i.e. without any metallic piece 
between the phase-differing superconductors, such modes 
are still present. On the other hand, when the width N,, 
increases we see a gradual smearing of the zero-energy 
quasiparticles and novel features appearing in the Majo- 
rana polarization. This is a consequence of the reduced 
proximity induced gap in wider strips which naturally 
reduces the localization of any mid-gap states. 
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FIG. 3. Components of the Majorana polarization Pio ob- 
tained for the region highlighted by the dashed frame in Fig. 2. 
The magnitude of the arrows shows |Pis;| and their direction 
shows Arg Pic. We note that the phase of the Majorana po- 
larization is only well defined up to a global shift. The shaded 
region is the metallic strip. 


B. Majorana polarization 


Majorana modes are quasiparticles with energy En = 0 
(we denote such doubly-degenerate eigenstate by n = ng) 
and which are eigenstates of the particle-hole transforma- 
tion operator. We analyze here another valuable source of 
information about these modes encoded in the Majorana 
polarization.?" 7945 This quantity is particularly useful 
for characterizing the Majorana modes of quasi two- 
dimensional topological superconductors?^?5:4! where its 
phase develops both longitudinal and transverse varia- 
tion. As we shall see, its texture brings an important 
message about the delocalized Majorana quasiparticles. 


Let us start by checking the contributions Pio from 
each spin o to the Majorana polarization in the narrow 
metallic strips, when the zero-energy quasiparticles are 
well localized near its ends. To be specific, we focus 
on the region marked by the dashed lines in the bot- 
tom panel in Fig. 2. Both constituents P;, are depicted 
in Fig. 3 on the lattice sites of the marked metallic re- 
gion. We clearly note that the directions of the arrows 
depicting P;+ are opposite to Pi, which is typical for 
both strictly one-dimensional topological superconduc- 
tors (see Fig. 2 in Ref. 40) and for higher dimensions 
as well.?^9?9:41 The magnitudes of the two components, 
shown by the length of the arrows, are however very dif- 
ferent. It is important to emphasize, that for the realiza- 
tion of a true MBS the local phase of Piy — Pi, must be 
constant. Such a constraint Piy = e? |P] = —e'^[P;| 
seems to be satisfied in our case only for narrow metallic 
strips. 


There are in fact two conditions on P; which are re- 
quired for Majorana modes. The first, which we have 
seen here, is that its phase must be constant. The phase 
of P; does not appear to be measurable in any simple 
way. What is measurable, as we will show in Sec. IV, 
is its magnitude |P;|. For a Majorana mode we require 
|P;| = pf, and this gives a measurable determination 
of Majorana modes. Upon increasing the width of the 
metallic strip the Majorana polarization gradually de- 
velops varying orientations, both along and across the 
sample. Emergence of the transverse gradient is very 
sensitive to the width Nw, as illustrated in Fig. 4. 


IV. POLARIZED ANDREEV SPECTROSCOPY 


Here we briefly discuss an empirical method based on 
spin polarized Andreev reflection spectroscopy, ? which 
could probe the absolute value of the Majorana polariza- 
tion. Let us consider a scanning tunnelling microscope 
(STM) tip approaching site i of our heterostructure. The 
influence of this external reservoir of itinerant electrons 
can be incorporated by augmenting the model Hamilto- 


nian Eq. (1) with the local term 


Hi; = ECT E Litip) Chg Cko = » (roche dic + Hc.) ] 
k,o k,o 
N. 


STM tip hybridization 

(11) 
where the chemical potential of the tip, uj = p + eV, 
can be varied by an applied bias voltage V. &y is the dis- 
persion of the tip electrons, and ^ the tunneling ampli- 
tude from the tip to the heterostructure and vice-versa. 
The quasiparticle spectrum at site i can be indirectly 
deduced from measurements of the charge transport in- 
duced by the voltage V applied between the STM tip and 
the sample. In the subgap regime, i.e. for e|V| € A, such 
a current would solely originate from Andreev scatter- 
ing processes. This mechanism relies on the conversion 
of electrons arriving from the STM tip into the Cooper 
pairs of the superconducting heterostructure, with holes 
being reflected back into the STM tip. 

We are interested here in probing the topological su- 
perconducting phase related to the p-wave pairing of 
equal spin electrons induced between neighboring sites, 
analogous to the situation in Kitaev's simple model.*° 
For this purpose let us assume a complete polarization 
of the tip, where only one spin component c can partici- 
pate in the charge transport. On a microscopic level we 
thus imagine that an itinerant electron of spin c arrives 
from the polarized STM tip at site à where it forms a 
triplet pair with an electron from the neighbouring site 
j, sending a hole of the same spin orientation into the 
tip. Within the Landauer formalism we can express the 
resulting spin-dependent charge current by^? 


I) = f dw T3(w) Uto eV) - flu ev) 
(12) 
where f(w) = [1--exp(w/kpT)]| | is the Fermi-Dirac 
distribution function. The main quantites of interest are 
the spatially-dependent transmission probabilities^? 


T£(w) =I} |F", (13) 


where Fj (w) = (dis; dla is the Fourier transform of 
the anomalous retarded Green's function in Nambu rep- 
resentation. For practical reasons, because we focus on 
a small charge transport window which is only a frac- 
tion of meV around the chemical potential u of our sam- 
ple, we have introduced a constant coupling strength, 
Iw = 2r X [k| E(w — ex). Formally it is equivalent to 
the wide band limit approximation. 

The anomalous Green's function Fg (w) has the ex- 
plicit form, for o =f, 


wn)" 
wtitly +E, }’ 
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where [Un and Uj, can be computed numerically. 
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FIG. 4. The Majorana polarization P; obtained for the heterostructure comprising Nw = 4 (top) and Nu = 10 (bottom) atomic 
rows in the metallic strip, marked by the shaded region. Numerical results are obtained for the same model parameters as in 
Fig. 2 but using Nz = 100, Ny = 20. The magnitude of the arrows shows |Pic| and their direction shows Arg Pic. We note 


that the phase is only well defined up to a global shift. 


Focusing on the zero-energy limit w — 0, dominated by 
the Andreev scatterings via the Majorana quasiparticle 
(Eno = 0), we can express the transmittance by 


c ps n no \* n no\*|2 
Ty(w~ 0) = STE jua (v2) + uso (vis) | » (15) 


Substituting Eq. (15) into the Andreev current formula 
(12) yields, at low temperatures, the following zero-bias 
differential conductance 

AV) den 


im ————— ~ 


2 
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This result (16) demonstrates that selective Andreev 
transport might probe the spin-dependent contribution 
(10) to the Majorana polarization. Strictly speaking, 
however, such tunneling processes occur on the links (in- 
volving the neighbouring sites i and j) rather than on 
individual local sites. For this reason the differential con- 
ductance (16) would measure the symmetrized Majorana 
polarization 

Pi» = Uig Ujg + Ujo Vio (17) 
over the neighboring sites i and j, instead of the strictly 
local definition (10). Since the diagonalization coeffi- 
cients are slowly varying in space, u$, © uj, Vj, © Vio» 
the symmetrized Pcij>,o and local Pio Majorana polar- 
izations should in practice be nearly identical. 

Let us finally recall that the complex vector uzy vy is 
typically perfectly opposite to Up UD. see Fig. 3. This 
implies, that the absolute value of the Majorana polar- 
ization (9) |P;| = |Pit| + [Pi], and the same holds for 
the symmetrized Majorana polarization. 
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FIG. 5. The spatial profile of |P;|/p? for the potential zero- 
energy Majorana quasiparticles appearing in the metallic strip 
consisting of 1, 2 and 5 rows of atomic chains, as indicated. 
We have used the same model parameters as in Fig. 2: A — 
0.25t, à = m, A = 0.5t, Bo = 0.1t, u = —3.75t. The difference 
between the well isolated MBS and the delocalized states in 
the wider strip is evident. 


In conclusion, by measuring the zero-bias differential 
dI?. : : 
conductance E of the spin-selective Andreev current 
flowing through the neighboring sites 7 and j one can 


evaluate the absolute value of the symmetrized Majorana 


polarization 
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As far as the spatial variation of the phase is concerned, 
its determination is, if possible, evidently more cumber- 
some. This problem is beyond the scope of the present 
study. 

Although the spatial variation of the phase is not cur- 
rently measurable we can compare the absolute value of 
|P;| with p;(0). For a MBS these must be the same, there- 
fore if we plot the ratio of these measurable quantities, 
|P;|/p;(0), it should be flat for a MBS. In Fig. 5 we show 
results for N,, € {1,2,5}, comparable to Fig. 2, which 
show that the MBS profile is indeed flat. For N,, — 5, 
when the two MBS at either end of the strip start to 
overlap, this quantity is no longer flat. Thus this can 
be used as an experimental determination of whether lo- 
calized states are actually MBS. It is worth noting that 
increasing the length Ns of the system would result in 
[P;|/p;(0) being flat even for Ny > 5. 


V. LOCALIZATION OF THE MAJORANA 
MODES 


In Sec. III we have shown that upon increasing the 
width N,, of a metallic strip the topological supercon- 
ducting state reveals (1) smearing of the Majorana quasi- 
particle, and (ii) development of the transverse gradient 
of the Majorana polarization. One may ask, however, 
whether there is any chance of localization of the Majo- 
rana modes. In this section we illustrate that such an 
effect could be observable due to local defects introduced 
in certain regions of the metallic strip. 

Topological superconductivity in 1-dimensional wires 
and atomic chains has been shown to be robust against 
weak disorder, 54 noise, inhomogeneous spin-orbit 
coupling,” reorientation of the magnetic field,* ^" and 
thermal fluctuations.°> ?! Sufficiently strong scatter- 
ing centers, however, could effectively break these 1- 
dimensional systems into separate segments, inducing ad- 
ditional pairs of the Majorana modes.9? Such a mecha- 
nism can be expected to be inefficient in 2-dimensional 
systems. To verify this conjecture for the quasi 2- 
dimensional heterostructure discussed in this paper we 
take into consideration a point-like electrostatic defect 
Himp = Vod! dic positioned at site 4g of the metallic 
region. We assume this scattering potential to be rather 
strong, Vy = 100t, as otherwise its influence would be 
less visible. 

Let us first assume the scattering potential to be placed 
in a central part of the metallic region (Fig. 6(a)). Nei- 
ther the zero-energy spectral function nor the Majorana 
polarization reveal any influence of such an electrostatic 
impurity on the existing Majorana quasiparticles, in con- 
trast to the properties of 1-dimensional topological super- 
conductors.9? This behaviour seems to be quite natural, 


because the Majorana modes are safely distant from the 
impurity. 

Contrary to this situation, let us next consider the scat- 
tering potential near the left side of the metallic strip 
(Fig. 6(b,c)). We selected the specific sites ij = 4 and 
ig = 8, in order to guarantee a considerable overlap of the 
local defect with the left hand Majorana mode. Under 
such circumstances the scattering potential has a sub- 
stantial influence both on the spectral function (left pan- 
els) and the Majorana polarization (right panels). This 
electrostatic impurity reduces the spatial extent of the 
Majorana quasiparticle on the left hand side, whereas 
the other Majorana quasiparticle is practically left in- 
tact. Such a tendency towards localization of the Ma- 
jorana modes has been recently predicted by Haim and 
Stern,’ when investigating the different role of weak ex- 
tended disorder. 

Our present study provides more detailed information 
concerning such disorder-induced-localization, indicating 
that: (i) disorder present in the internal segments of the 
metallic strip would naturally be rather ineffective for 
the Majorana quasiparticles, whereas (ii) disorder intro- 
duced to the regions of already existing Majorana quasi- 
particles substantially reduces their spatial extent. The 
considerations discussed in Sec. IV suggest that empiri- 
cal detection of this subtle phenomenon could be feasi- 
ble. A tendency towards the localization of the Majorana 
quasiparticles could be observed in the maps of differen- 
tial conductance for the spin polarized Andreev current 
induced via the metallic region in the presence of the 
intentionally deposited local defects. Such an electro- 
static scattering potential could be created by applying 
gate potentials, whereas a magnetic potential, leading to 
similar effects, can be obtained by locally perturbing the 
Zeeman field. 


VI. SUMMARY 


We have theoretically studied the properties of the Ma- 
jorana quasiparticles emerging in a narrow metallic strip 
sandwiched between two s-wave superconductors in a 
Josephson-junction geometry. The topological supercon- 
ducting phase has been recently reported for such metal- 
lic strips by the groups in Copenhagen?? and Harvard?’ 
with a length-to-width ratio ranging from 20 to 100, re- 
spectively. Using the Bogoliubov de-Gennes treatment 
we have investigated the role of the finite metallic strip 
width, exploring its influence on: (a) spatial profiles of 
the zero-energy quasiparticles and (b) topography of the 
Majorana polarization that probes the particle-hole over- 
lap of the zero-energy quasiparticles. Furthermore, we 
have proposed a feasible method for detecting the mag- 
nitude of such a Majorana polarization by measuring the 
differential conductance in spin-polarized Andreev reflec- 
tion spectroscopy. 

We have also analyzed the influence of strong (point- 
like) electrostatic defects on the Majorana modes. We 


FIG. 6. The spatial density profile of the zero-energy Majorana quasiparticles (left panels) and their polarizations (right panels) 


obtained in the proximitized metallic strip, consisting of N, 


= 5 rows of atoms, in presence of a point-like electrostatic defect 


with Vo = 100t at site: (a) io = 40a, (b) io = 8a, and (c) io = 4a and Ny = 20. The magnitude of the arrows in the right hand 
side plots show |P;;| and their direction show Arg Pic. We note that the phase is only well defined up to a global shift. 


have revealed that such a local scattering potential can 
affect the localization length of the Majorana quasiparti- 
cles if deposited near the ends of the metallic strip. Under 
such circumstances the neighboring Majorana mode sub- 
stantially reduces its spatial extent, which can be com- 
pared with what has been predicted in Ref.,°° whereas 
the opposite-end Majorana mode remains practically in- 
tact. Similar coexistences of the localized and delocalized 
Majorana quasiparticles have been previously observed 
by scanning tunneling spectroscopy using a disordered 
monolayer of superconducting Pb coupled to underly- 
ing Co — Si magnetic islands.?? We hope that proxim- 
itized metallic strips would be a convenient platform not 
only for the realization of topological superconductivity, 
tunable by the magnetic field and Josephson phase, but 
could also allow for manipulating the Majorana quasi- 
particle length-scale. 
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FIG. 7. The topological superconducting phase (dark area) 
with respect to the magnetic field Bo and the Josephson 
phase $ obtained by the criterion concerning the global Ma- 
jorana polarization.?^?5 Numerical computations were done 
for Nw = 2 and the same model parameters as in Fig. 2. 


Appendix: Influence of the Josepshon phase 


In the main part of this manuscript we have analyzed 
the spectroscopic properties of the Majorana quasipar- 
ticles, focusing on the particular case ¢ = 7, that is 
optimal for occurrence of the topological superconduct- 
ing state. Similar qualitative features would be also ob- 
servable for other values of the Josephson phase ¢ Æ 7, 
provided that model parameters (such as, e.g., the mag- 
netic field) are appropriately tuned.?9?"7 A possible crite- 
rion (proposed by one of us???) for determination of the 
topological diagram relies on the Majorana polarization 
for MBS to be non-vanishing when summed over a por- 
tion of the investigated system where the bound states 
should reside, in fact it should be equal to the total den- 
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